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SIMULATION AND ANALYSIS OF A GEOPOTENTIAL RESEARCH MISSION 


Abstract 


The Geopotential Research Mission (GRM) is a proposed mission that uses two 
satellites in nearly identical orbits to determine the Earth's geopotential field using 
measurements of relative range rate between the two satellites. In the current proposal, 
the satellites will be in polar orbits, at an altitude of 160 kilometers, but separated by 
several hundred kilometers. A drag compensation system will be used to eliminate the 
nonconservative forces that affect the purely gravitational trajectories of the satellites. 

This study has investigated methods for the determination of the initial 
conditions for the two satellites that will satisfy the mission requirements. For certain 
gravitational recovery techniques, the satellites must remain close to a specified 
separation distance and their groundtracks must repeat after a specified interval of time. 
Since the objective of the GRM mission is to improve the gravity model, any pre- 
mission orbit predicted using existing gravity models will be in error. A technique has 
been developed to eliminate the drift between the two satellites caused by gravitational 
modeling errors and return them to repeating groundtracks. This technique is 
independent of the geopotential field and other perturbations that might have been 
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neglected in the pre-mission model. 


The concept of "frozen orbits" was also investigated for possible GRM orbits. 
A frozen orbit restricts the secular motion of the argument of perigee and removes the 
long period changes of the eccentricity. This characteristic of the frozen orbit 
minimizes the altitude variations over given latitudes on the Earth. Frozen orbits also 
have the further advantage of more easily maintaining a repeating ground track. 

The effects of temporal perturbations on the relative range-rate signal are also 
investigated. At the proposed altitude of 160 km, the range-rate signal produced by 
perturbations other than the static geopotential field are dominated by the luni-solar 
effect. This study demonstrates that the combined effects of all the temporal 
perturbations does not prevent the orbit from being frozen, nor do they prevent the 
satellites from obtaining a repeating groundtrack to within the specified closure 
distance. 


IV 



TABLE OF CONTENTS 


Abstract iii 

1. Methods for the Recovery of the Geopotential Field 

1.1 Introduction 1 

1.2 Review of Gravity Models 2 

1.3 Ground Based Satellite Tracking 4 

1.4 Satellite Altimetry 5 

1.5 Gravity Gradiometers 7 

1.6 High-Low Satellite Pairs 9 

1.7 Low-Low Satellite Pairs 11 

1.8 Introduction to Topics 12 

2. Concept of GRM 

2.1 Introduction 14 

2.2 Prior work on GRM 16 

2.3 Description of Mission 16 

2.4 Mission Requirements and Goals 19 

2.5 Specifications for Mission Simulation 21 


v 


3. Frozen Orbits 

3.1 Background on Frozen Orbits 24 

3.2 Definition of Frozen Orbits 26 

3.3 Characteristics of Frozen Orbits 28 

3.4 Addition of Short Period Terms 31 

3.5 Non-Frozen Orbit 33 

3.6 Summary 35 

4. Adjustment of Initial Conditions 

4.1 Introduction 51 

4.2 Determination of a Principal Set of Initial Conditions 54 

4.3 Use 32 days to Determine Initial Conditions 56 

4.4 Non-Polar Adjustments 62 

4.5 Sensitivity Study 62 

4.6 Range of Linear Reliability 67 

4.7 Correction of Velocity 69 

4.8 Earlier Determination Orbital Adjustments 71 

4.9 Summary 77 

5. Analysis of GRM simulation 

5.1 Introduction 97 

5.2 Description of Simulation 98 

vi 


5.3 Investigation of General Behavior 103 

5.4 Reduction of Residuals 109 

5.5 Summary 115 

6. Effect of Temporal Perturbation on Relative Motion 

6.1 Introduction 159 

6.2 Precession, Nutation and Polar Motion 160 

6.3 Solid Earth Tides 163 

6.4 Ocean Tides 164 

6.5 Planetary Effects 165 

6.6 Luni-Solar 166 

6.7 Relativistic Effects 168 

6.8 Effects of the Perturbation on the Initial Conditions 169 

6.9 Effect of the Perturbation on the Frozen Orbit 170 

6.10 Summary 171 

7. Conclusion/Summary 

7.1 Summary 192 

7.2 Additional Research 194 

Appendix A 196 

Appendix B 199 

Bibliography 201 

vii 


CHAPTER 1 


METHODS FOR THE RECOVERY OF THE GEOPOTENTIAL FIELD 
1.1 Introduction 

To improve the present knowledge of the Earth's gravity field, a mission 
dedicated to the recovery of the geopotential field has been suggested. In 1978, the 
National Academy of Science Committee on Geodesy suggested that two, low altitude, 
polar satellites be studied as a method for the recovery of the Earth's gravity field 
[National Academy of Science, 1979]. The primary goal of this mission, currently 
known as the Geopotential Research Mission (GRM), would to produce a high 
resolution global gravity field model of the Earth. 

In the past, modem geopotential models of the Earth were tailored for a 
particular satellite, thereby producing accurate results when tested on that particular 
satellite, but less accurate results when applied to other satellites with different orbit 
characteristics [Rapp, 1981]. There is a need for a global gravity model, not tuned for 
a particular satellite, but one that would be useful for both terrestrial and satellite 
applications. The rationale for the improvement in gravity field models will be 
explained in this chapter. 

The Earth's gravity field is nearly complete between ±70° latitude with a 
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resolution of 1° [Lerch, Putney, Wagner and Klosko, 1981]. However, the best 
resolution is mainly in the ocean areas, while much of the solid Earth area is 
inadequately determined. Also, the shorter wavelength harmonics of the geopotential 
field need to be discerned in order to observe the more detailed structure of the Earth's 
surface and the mechanics of its motion. An accurate knowledge of the gravity field 
will provide information on the internal structure of the Earth, as well as the surface 
features. With improvements in the gravity model comes an improvement of the 
Earth's geoid. The geoid is the mean sea level, or the shape the Earth would be if there 
were no land masses and no tidal effects; it is a surface of constant gravitational 
potential [ Stewart , Lu and Lefebvre, 1986], An accurate geoid provides clear 
information on the general shape of the Earth. 

Though more progress through better use of data already acquired can be made, 
new techniques and refinements of old ones need to be devised for the improvement of 
the global geopotential field. Since current methods for measuring the Earth's gravity 
field are incapable of providing a high resolution, homogeneous, global geopotential 
model [ESA, 1987], the proposed satellite mission dedicated to the recovery of the 
Earth’s gravity field will use an alternate method to those employed in the past. 
Various methods have been proposed for this mission. The candidate methods include 
satellite-borne gravity gradiometers, and satellite-to-satellite range-rate measurements 
for high-low and low-low satellite configurations. The advantages and disadvantages 
of each of these new methods along with a description of current techniques are 
discussed in the following sections. 
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1.2 Review of Gravity Field Models 

Prior to the first Earth satellite launched in 1957, the determination of the 
Earth's gravity field was restricted to surface gravimetry measurements [Lambeck and 
Coleman , 1983]. With the use of satellites, more measurement data from various other 
techniques could be applied to the knowledge of the Earth's geopotential field. Three 
representative gravity models that are used in this study and that will be discussed in 
this section are the Goddard Earth Model GEM10B, and the models provided by Ohio 
State University, referred to as the OSU322 and the OSU86F fields. 

For the past fifteen years, Goddard Earth Models have been under development 
[Lambeck and Coleman, 1983]. The odd numbered gravity models computed by the 
Goddard Space Flight Center, such as GEM 9, are based strictly upon satellite tracking 
data. The even-numbered models differ from the odd-numbered models by the 
inclusion of surface gravimetry data as well [Lerch, et al., 1981]. The accuracy of 
these fields was assessed by Lerch, et al. [1981]. 

The Geodynamic Experimental Ocean Satellite, GEOS-3, the first unmanned 
satellite to provide altimetry data, was an important source of information for the most 
recent Goddard Earth Models. In 1977, the altimetry data from GEOS-3 was used to 
improve the Goddard Earth Model GEM 10 from a field complete through degree and 
order 22 (22 x 22) to the more current GEM 1 OB, a 36 x 36 geopotential field [Lerch 
and Wagner, 1981]. 


Ohio State University (OSU) has been developing gravity models for the past 
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twenty years [Lambeck and Coleman, 1983]. In the creation of the latest OSU fields, 
Rapp [1981] attempted not to tune the fields for a particular satellite, but instead 
produce fields that were more representative of many satellites. The OSU322 field is 
180 x 180 with some modifications that include terms to degree 300 [Rapp, 1981], and 
OSU86F is a full 360 x 360 field [Rapp and Cruz, 1986]. Though these fields have 
been expanded to over degree 180, the high degree and order harmonic coefficients 
were determined mainly from surface gravimetry data [Rapp, 1981]. 

To obtain the finer structure of the Earth's gravity field, surface gravimetry data 
must be included along with the satellite tracking observations. However, 100% errors 
in the terms over degree 150 have been found. In addition, surface data does not cover 
the entire globe, and because the actual acquisition of the data tends to vary, the data 
from surface gravimetry measurements is not uniformly distributed [ESA, 1987]. 

The accuracy to which a satellite's orbit can be determined is sometimes used as 
a criterion for evaluating certain geopotential models. It is desired, however, to have a 
gravity field model that is a physical representation of the Earth, not one that produces 
accurate motion for an individual satellite or even a group of satellites in similar orbits 
[Lambeck and Coleman, 1983]. 

1 . 3 Ground Based Satellite T racking 

Until the introduction of satellite altimetry, models of the Earth's gravity field 
were derived almost exclusively from ground based satellite tracking data and from 
surface gravimetry data [Lerch, et ai, 1981]. Satellite tracking data provide the longer 
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wavelength features of the gravity field, whereas the gravimetry data provide the finer 
details [Smith, Lerch, March, Wagner, Kolenkiewicz and Khan, 1976]. 

A satellite’s sensitivity to certain geopotential coefficients is dependent upon its 
orbit; consequently, to properly determine the gravity field, data from several satellites 
that are at different altitudes and varying inclinations should be obtained [Lambeck and 
Coleman , 1983]. In order to track the satellite's orbit, ground based tracking stations 
must be available, and collectively must be able to track the satellite long enough to 
enable determination of an accurate orbit. From satellite tracking observations only, the 
gravity field can be determined up to degree and order 20 plus some additional zonal 
harmonics [ESA, 1987], although the accuracy of the coefficients depends also on the 
precision of the observations. 

Over the years, the global satellite tracking networks have become more 
extensive. Tracking networks have been established by the U.S. Navy (known as 
TRANET), by the Smithsonian Astrophysical Observatory (SAO), and by Goddard 
Space Flight Center (GSFC) [Torge, 1980]. These networks have more than a total of 
200 stations around the globe, and are capable of measuring a satellite's state using 
laser, optical, electronic, and doppler observations. 

Due to physical and political constraints, the ground based stations are not 
placed uniformly around the Earth and, therefore, some geographic areas have very 
sparse tracking data [Ar gender o and Lowery, 1978]. The gravity models are more 
accurate in areas where stations are located and large amounts of data can be 
accumulated. The lack of global coverage and uniformity in tracking is a significant 
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disadvantage of ground based tracking for the determination of a high accuracy and 
resolution gravity field. The limitations indicate that additional measurement techniques 
for the recovery of the global geopotential field must be employed. 

1.4 Satellite Altimetry 

The satellite altimeter measures the distance between the satellite and the mean 
ocean surface directly beneath it. The altimeter signal is transmitted from the satellite to 
the ocean surface and received back at the satellite. The time required for the signal to 
be sent and received again is a direct measurement of the satellite altitude above the 
ocean surface. Because the geoid is a "global reference surface for height", an accurate 
geoid is essential for some applications of altimetry measurements [Torge, 1980]. With 
precise altimetry data, more detailed information on ocean circulation and the currents' 
motion can be obtained. By comparing the ocean surface as measured by the altimeter 
to the geoid, ocean currents can be observed. 

Until the first altimeter was flown on Skylab on January 31, 1974, the Earth's 
gravity models were based mainly on ground based satellite tracking data and on 
surface gravimeter measurements. Current gravity models have been determined using 
not only satellite tracking and surface gravimeter data, but also altimeter data from one 
or more satellites such as GEOS -3 and SEAS AT. The use of altimeter data from these 
two particular satellites has further improved the ocean surface accuracy level to the two 
meter range and better in some areas. 


In addition to biases caused by drifts in the satellite's clock, random errors in 



7 


the altimeter measurements occur from thermal noise and from instrument limitations. 
The atmosphere and ionosphere, as well as the ocean surface will also introduce errors 
depending on the signal frequency and the roughness of the sea. Theoretically, all the 
errors in single frequency radar measurements can be reduced to approximately the 8 
cm level [Greene, 1971]. SEASAT demonstrated a precision level of 10 cm or better 
[Lerch, Marsh, Klosko and Williamson , 1982]. Current technology proposals include 
the use a dual frequency radar altimeter signal to reduce measurement errors by 
providing a direct measure of the ionospheric effects. TOPEX will use this technique, 
and is expected to have the altimeter measurement errors to within 2 cm [Stewart, et al., 
1986]. To achieve this level of accuracy for altimeters, the knowledge of the 
geopotential field must be improved significantly [Lerch, et al., 1982]. 

Altimetry is mainly used for the recovery of the shorter wavelength features of 
the gravity field, but it is not a very useful technique for the recovery of the longer 
wavelengths [ESA, 1987]. Large discrepancies in the geopotential field remain in both 
the land areas and in the ocean areas above ±72° latitude where GEOS-3 and SEASAT 
were unable to cover [ESA, 1987]. Altimetry data are useful for improving the gravity 
field knowledge of the oceanic areas, but not the continental regions of the globe. Even 
with the combination of satellite tracking data, surface gravimetry data, and altimetry 
data the resolution of the gravity field cannot be obtained to the accuracy levels that are 
desired, and therefore other methods must be applied [ESA, 1987]. 

1.5 Gravity Gradiometers 


Space-borne gradiometers have been under development since 1970, although 
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to date, no gradiometers are known to have flown in a satellite mission. However, 
gradiometers have been used to measure the gravity field from instruments carried on 
airplanes, as well as on the ground [Argentiero and Lowrey, 1978]. The most recent 
development in moving based gravity gradiometers is the Gravity Gradiometer Survey 
System (GGSS). The GGSS, consisting of three pairs of mutually orthogonal 
accelerometers, is installed in a mobile van that can measure the gravity gradient by 
either traveling on the ground or by being carried in a C- 130 airplane [Jekeli, 1987], 

In the presence of a gravitational field, the gravity gradiometer measures the 
differences in the forces sensed by two or more accelerometers. The amplitude of the 
output signal from the sensors supplies the second derivatives of the gravitational 
potential from which the gravity values can be estimated [Forward, 1972]. The 
uncertainties in the measurements have been improved in the last few years from less 
than one Eotvos Unit (E.U.) in 1980 [Torge ], to between 10' 3 E.U. and 10" 5 E.U. 
[ESA, 1987], where one E.U. is equivalent to 10' 6 m s" 2 /km [Torge, 1980]. 

The long wavelength features of the gravity field are relatively well known, 
therefore, the recovery of the short wavelength terms in the gravity field is of principal 
importance for obtaining a comprehensive gravity model. These terms have little 
contribution to the overall gravity field, but can contribute significantly to the gravity 
field in the neighborhood of an anomaly. The gradiometer's advantage is that it senses 
the gradient of the gravity field of the immediate area, enabling it to resolve the short 
wavelength terms [Forward, 1972]. 


Feasibility studies have been performed on the usefulness of gradiometers in the 
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recovery of the gravity field. For a satellite at an altitude of 200 km, a 1° x 1° gravity 
anomaly of 1 to 2 mgal will require the gradiometer to have a measurement accuracy of 
under 10" 2 E.U. [ Reigber , Keller, Kunkel,and Lutz, 1986]. The results from studies 
performed by European agencies indicate that gravity fields determined with 
gradiometers are comparable in accuracy to those obtainable from alternate satellite 
measurement types. The proposed Gravity Gradiometer Mission (GGM) anticipates an 
accuracy in the measurement signal of 10^ E.U. This accuracy level would produce a 
50 km resolution of the gravity field [Paik, 1985]. Unlike current altimetry data, 
however, the gravity gradiometers can provide global data for the gravity field, without 
restriction to specific geographic areas. Methods for the recovery of the gravity field 
from satellite gradiometry are at a very early stage. A preliminary study of the 
determination of the gravity field using full tensor gradiometry was described by 
Colombo [1987]. The primary disadvantage of the satellite gradiometry technique is 
that it is not a proven or well-tested technique for satellite applications. 

1.6 High-Low Satellite Pairs 

The "high-low" satellite configuration is a proposed technique for meeting the 
dedicated gravity mission requirements. It consists of one high altitude satellite and 
another satellite in a much lower altitude, polar orbit. The low satellite must be 
sufficiently low to sense short wavelength gravity anomalies, while the high satellite 
must be sufficiently high to be insensitive to these short wavelength features of the 
gravity field, as well as atmospheric effects [S/ry, 1973]. This concept was originally 
proposed to recover long wavelength effects, but it can be used for the detection of the 
short wavelength terms of the gravity field as well [Agentiero andLowrey, 1978]. 
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In 1972, the high-low configuration of the ATS-F/Nimbus-E was suggested as 
an experiment to recover the short wavelength features of the geopotential field [Von 
Bun, 1972]. ATS-F, a geostationary satellite, and Nimbus-E, a weather satellite at an 
altitude of 1000 km, demonstrated that the range-rate data from satellite-to-satellite 
tracking could successfully be used to detect gravity anomalies of short wavelengths 
even though the Nimbus satellite had a fairly high altitude for this type of recovery 
technique. With the encouraging results from ATS/Nimbus, a high-low mission 
dedicated to recovering the gravity field was proposed. 

The Geopause/Gravsat mission was first proposed in 1973. The high altitude. 
Geopause satellite would have been a 14 and 24 hour, polar orbit. The low altitude, 
Gravsat satellite, was proposed to have an approximate height of 300 km, in a polar 
orbit as well. Polar orbits were selected in order to have complete global coverage of 
the gravity field, and the 300 km altitude was regarded as low enough to sense the short 
wavelength gravity perturbations. To maintain the satellite at the desired altitude, the 
low satellite would have to be equipped with a drag compensation device to adjust for 
atmospheric effects that would cause the orbit to deteriorate. This compensation device 
would allow the satellite to follow a purely gravitational path. 

With both the high and low altitude satellites in the same orbit plane, the 
Geopause satellite could observe the along-track and radial components of the Gravsat 
orbit but, it would be less sensitive to the cross-track component By adjusting the two 
orbit planes to have an angle of 30° to 45° between them, it was shown by Siry [1973] 
that the cross track component could be observed by the higher orbit satellite. 



In 1974, Koch and Argentiero [1974] performed a simulation of the 
Geopause/Gravsat mission. Their main interest was to improve the lower degree and 
order coefficients of the gravity field. They expected one or two orders of magnitude 
of improvement of the terms up to degree 8 and order 6. However, Estes and 
Lancaster [1976a] demonstrated that the same mission could be used to recover the fine 
structures of the gravity field, and by lowering the altitude of the Gravsat satellite, the 
short wavelength terms could be more easily detected. 

1.7 Low-Low Satellite Pairs 

As early as 1969, it has been proposed that a satellite pair, in identical polar 
orbits, could effectively recover the Earth's gravity field [Wolff, 1969]. The basic idea 
of this "low-low" configuration is that the satellites would maintain a low, polar orbit 
that should completely cover the entire Earth in approximately one month's time. These 
satellites would be separated by a specified distance which would depend upon the 
selected altitude and the level of accuracy desired. The measurement signal would be 
the relative motion, either range or range-rate, between the two satellites. 

In addition to the high-low mission simulation previously discussed, Estes and 
Lancaster [1976b] also performed a simulation for the low-low mission. The altitude 
of the satellite pair was specified to be 250 km. This study showed that the low-low 
configuration would sense the local gravity field (much like the gravity gradiometers) 
and would be less sensitive to gravity effects that are further away than the conventional 
satellite measurement signal, since anomalies would affect both satellites in a similar 


manner. Their overall results, however, were not very promising. There was 
difficulty in the measurement of both the radial and the cross track components of the 
relative velocity between the two satellites. They also found a significant problem with 
differentiating between the frequencies generated by the gravity field, i.e. aliasing, and 
they concluded that the high-low configuration was more favorable. 

The most significant advantage of the low-low configuration over the high-low 
configuration is the length of mission time required to completely survey the Earth. 
The low-low mission would need only four weeks to completely survey the Earth, 
whereas, the high-low mission would require four months to obtain the same coverage 
of the Earth [ESA, 1987]. The shorter time needed for the low-low mission allows for 
the same mission to be performed repeatedly. Variations of the mission can be 
performed by changing the separation distance in order to cause the measurement signal 
to be more sensitive to certain wavelengths. This ability to repeat the mission will 
insure the data integrity and will help to eliminate the aliasing problem described by 
Estes and Lancaster [1976b]. Another important advantage of the low-low mission is 
that the recovery of the geopotential coefficients are expected to be one order of 
magnitude better than could be achieved with the high-low mission [Willis and Smith, 
1980]. 


The satellite-to-satellite methods have fewer measurement error sources than the 
ground based tracking methods. Both ground station locations and the atmosphere are 
major sources of errors in satellite measurement techniques, however, the low-low 
satellites always remain in sight of each other and would have continuous global 
coverage, a major advantage over ground based measurements. The low-low 


configuration was selected for study in the Geopotential Research Mission because of 
its greater sensitivity to the gravity anomalies [ESA, 1987], and because of the short 
mission time it requires. A discussion of a proposed concept for this mission will be 
described in more detail in the following chapter. 

1.8 Introduction to Topics 

The objective of this dissertation is to investigate the aspects of the Geopotential 
Research Mission using the low-low configuration that pertain to mission planning and 
orbital operations. Requirements of the gravity field recovery phase of the mission are 
discussed. A simulation of the GRM satellite pairs was used to study the feasibility of 
the specified requirements. 

In Chapter 2, a discussion of the GRM satellite pair is presented along with the 
mission goals and expectations for the geopotential field recovery. In Chapter 3, the 
characteristics of frozen orbits are studied and are compared to non-frozen orbits. 
Methods for determining the initial conditions for the mission planning and the 
operational phase are described in Chapter 4. A simulation of the GRM satellite pair is 
studied and an identification of significant resonant terms are made in Chapter 5. In 
Chapter 6, the effects of temporal perturbations on the relative range-rate between the 
satellites are investigated. 


CHAPTER 2 


CONCEPT OF THE GEOPOTENTIAL RESEARCH MISSION 
2.1 Introduction 

The Geopotential Research Mission (GRM) is part of a program designed to 
recover the higher order geopotential coefficients by using the relative range-rate 
changes between two satellites in nearly identical orbits. The two satellites should 
maintain a nearly constant separation distance, typically several hundred kilometers. 
The orbits will be low altitude, polar orbits and the groundtracks should repeat after a 
specified number of days (Figure 2.1). The altitudes of the satellites depend upon the 
number of days required for the groundtrack to repeat, and on the number of exact 
revolutions of the satellites required in that same amount of time. To date, the altitude, 
the frequency of the groundtrack repeat, and the distance separating the satellites have 
not been finalized for the GRM, but a nominal mission is provided in this discussion. 
The mission is designed to have an operational lifetime of a minimum of six months 
[Keating, Taylor, Kahn and Lerch, 1986]. 

Spatial variations in the Earth's geopotential field can be measured by the 
changes in the relative range-rate between the two satellites. As the leading satellite 
approaches a gravity anomaly, it will respond with a change in its absolute velocity. 
This change is sensed through the doppler shift in the radar signal by the second, 
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trailing satellite, which has not yet been affected by the anomaly [NASA, 1984], 
Relative velocity is used as the measurement type because it is more sensitive to the 
variations in the gravity field than the relative range measurement [ESA, 1987]. 

Along with gathering information on the gravity field, these low altitude 
satellites are also expected to gather information on the Earth's magnetic field. The 
magnetic field would be measured by magnetometers; both scalar and vector 
magnetometers are to be carried by the leading spacecraft. The magnetometers will be 
placed on a boom that will shield them from the magnetic effects created by the satellite 
itself [Keating, etal., 1986]. 

The primary advantage of satellite-to-satellite tracking is the availability of a 
nearly continuous measurement signal. Due to geophysical and political constraints, 
the data from previous techniques used to recover the Earth's geopotential field are 
unevenly distributed. The limits of useful information from current techniques have 
essentially been achieved and new techniques must be devised to obtain a more accurate 
gravity model. The low-low configuration the chosen configuration for reasons 
discussed in Chapter 1 [Smith, Langel and Keating, 1982]. 

The ability to recover the geopotential field depends upon altitude, separation 
distance, data type, data rate, data noise, model errors and a priori values of the 
geopotential coefficients [Estes and Lancaster, 1976a]. The limitations and 
requirements on these parameters are discussed in this chapter, except for the a priori 
values of the geopotential, which will be discussed in Chapter 5. The research 
described in this dissertation will concentrate on the gravitational aspects of this 
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mission, but not upon the actual recovery of the global geopotential field. 

2.2 Prior Work on the Geopotential Research Mission 

A study presented by Colombo [1984] provided a set of initial conditions which 
produced an orbit that met many of the requirements designated for the mission. The 
trajectory from these initial conditions were compared to the trajectory from a new set 
of initial conditions derived for this study. This dissertation concentrates, in part, on 
determination of the best set of initial conditions that meet the requirements specified for 
this mission. The initial conditions are then used in a simulation for 32 days of mission 
lifetime. 

The GRM study is coordinated by the NASA Goddard Space Flight Center. 
Most of the previous work on GRM has concentrated on the recovery of the 
geopotential field. Wagner and Goad [1982], Kaula [1983] and Colombo [1984] have 
each devised a method to recover a high resolution global gravity field of the Earth. 
Results from the simulation of the GRM satellites presented in this study may be used 
to test these recovery techniques. 

2.3 Description of the Mission 

Originally planned to be launched from the shuttle, the orbit insertion plan for 
the GRM satellites may change in light of the recent problems with the shuttle. 
Regardless of the transportation mode into space, both satellites will be inserted into 
polar orbits at an altitude of 275 km with an initial separation distance of 50 km. The 



satellites will then undergo a series of partial descents and checkouts for four days until 
they reach the specified altitude of 160 km and a separation distance of 300 km. At that 
time, :nere will be seven days of final testing before the satellites begin their operations 
[Keating, et al., 1986]. Therefore, there will be only one week to place the two 
spacecraft into their proper orbits before the initiation of the operational phase of the 
mission. 

Since the recovery of the higher harmonics in the Earth geopotential field is a 
primary goal of this mission, the altitudes of the satellites should be as low as possible 
to enhance to sensitivity of the gravity signal. Atmospheric drag problems are avoided 
by using a Disturbance Compensation System (DISCOS) to cancel the effects of drag. 
First used on the TRIAD satellite in 1972, DISCOS contains sensors that activate 
thrusters in order to maintain a proof mass inside a cavity within a specified zone. This 
system not only shields the satellites from atmospheric perturbations, but also negates 
the effects from other nongravitational perturbations, such as solar radiation pressure, 
that may affect satellite motion. Enough fuel must be carried to compensate for the drag 
forces so that the satellites can maintain their orbits for at least six operational months 
[NASA, 1984]. The fuel requirements are an important factor in determining the 
limiting altitude of the satellites. 

The average altitude of the satellites is proposed to be approximately 160 
kilometers above the reference ellipsoid. If the altitude is lower than 160 km, then the 
satellites will have to correct more frequently for atmospheric effects, thus, requiring 
more fuel. A difference of only ten kilometers results in as much as a 35% increase in 
fuel consumption [Ray, Jenkins, DeBra and J unkins, 1985]. The altitude cannot 
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remain precisely constant at 160 km because of short period fluctuations in the 
semimajor axis which are caused by the perturbing forces. However, the average value 
of the altitude should remain near 160 km. The variation in altitude resulting from the 
gravitational forces will affect the DISCOS system, which will need to compensate for 
the short period changes in the drag forces associated with altitude variations. 

The semimajor axis is selected to produce a repeating groundtrack after a 
specified number of days, for example 30, 60, 90 or 180 days, depending on the 
desired equatorial spacing [Keating, et al., 1986]. The frequency of the groundtrack 
repeat chosen will maintain the satellites at a mean height above the reference ellipsiod 
of approximately 160 kilometers. In addition, the longitude of the ascending node is 
equal to 90° which places the orbit in the inertial Y-Z plane, thereby minimizing the 
luni-solar effects [Estes and Lancaster, 1976b]. The mean argument of perigee is equal 
to 90°, and since the orbit is polar, the inclination is equal to 90°. The eccentricity is 
selected so that it will have no long periodic effects and so there will be no long 
periodic or secular effects in argument of perigee. Such an orbit is referred to as a 
"frozen" orbit [Cook, 1966]. 

With a frozen orbit, the mean orbital ellipse does not change its shape or 
orientation, except for precession of the longitude of the ascending node, which will be 
close to zero since the orbits are polar. The orbits require constant mean orbital 
elements in order to allow the altitude above a particular point over the Earth to remain 
essentially constant with each satellite pass. This is a necessary condition for certain 
recovery techniques that use the Fourier harmonics in the determination of the 
disturbing function. Frequencies of the perturbations due to the geopotential will be 
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constant if the mean argument of perigee is also a constant [Wagner and Goad, 1982]. 
In addition, the small value for the eccentricity needed for the frozen orbit limits the 
variation in the semimajor axis. This fact restricts rapid changes in the altitude and will 
reduce the tracking errors for the individual satellites. 

Relative range-rate has been selected as the measurement signal because it is 
more sensitive than relative range to the changes in the local gravity field. The lead 
satellite reacts to a gravity anomaly before the second satellite does. This creates a 
change in the relative range-rate which is measured by the tracking systems onboard 
both satellites. When the trailing satellite reaches the gravity anomaly, there is another 
change in relative range-rate. 

Two doppler signals are continuously sent and received between the two 
satellites at the 42 and 91 GHz frequencies. These signals provide a data type known 
as one way integrated doppler. The two frequencies are required to correct for 
ionospheric refraction [NASA, 1984]. The signal broadcasted by the individual 
satellite is independent of the signal it receives. The combination of the signals between 
the satellites produces the final measurement, the relative range-rate. The shift in the 
doppler frequency determines the change in the relative range-rate, which in turn, is a 
measurement of the strength of the gravity anomaly [Keating, et al., 1986]. 

2.4 Mission Requirements and Goals 

The goal for the GRM is to obtain and improve mathematical models for the fine 
structure of the geopotential and the magnetic field. To be of geophysical interest, the 
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accuracy level for the gravity field improvement is required to be 2.5 mgal. This 
improvement should enable the geoid to be measured to within 10 cm while the 
magnetometers should measure the magnetic field to 1 nT (nanotesla). Both the gravity 
field and the magnetic field should be determined to a resolution of 100 km. The geoid 
is currently known from 20 to 50 cm with a resolution of 100 to 200 km, depending on 
the geographical location and on the gravity model [Smith, Langel and Keating, 1982]. 

The strength of the gravity signal decreases as the altitude increases, so an 
altitude increase would either limit the recovery of the higher degree and order gravity 
harmonics, or would require more precision in the relative range-rate measurement. At 
a 200 km altitude, the amplitude of the gravity signal decreases by a factor of 1.5 from 
one at a 160 km altitude; at a 250 km altitude, it decreases by a factor of three. The 
altitude must be low enough to enable resolution of the 2.5 mgal signal from a gravity 
anomaly that is 1° x 1° with the relative range-rate sensor precision of ±1 |im/s. The 
Goddard Earth Model (GEM 10C) has an expected accuracy of 20 mgal over a 1° x 1° 
area [Lerch, et ai, 1981]. The 160 km altitude produces stronger gravity signals than 
the higher altitudes and detects higher gravity harmonics at a given signal error. If the 
altitudes of the satellites were higher, then the higher degree harmonics would be too 
weak to be detected within the accuracy limits selected [Kahn and Felsentreger, 1982]. 
The gain in gravity signal strength is linear with the decrease in altitude, but the effect 
of drag increases exponentially with the same lowering of height, thus increasing the 
fuel requirements [Lowrey, 1975]. The higher the altitude, the less the capability to 
differentiate between the gravity anomalies, i.e., a loss of resolution occurs. As the 
altitude is reduced, the correlation between coefficients is decreased and there will be a 
gain in statistical independence between the harmonics [Estes and Lancaster, 1976b]. 
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At this point, it should be clear that the separation distance between the satellites 
is relevant to the detection of the gravity anomalies. Features on the Earth's surface that 
cause a change in the relative range-rate must be smaller than the separation distance 
between the satellites in order for the individual satellites to react to the anomaly at 
separate time intervals. The nominal separation distance of 300 km was selected in 
conjunction with the altitude. To obtain a variation in the measurement data generated 
by the gravity field, this value could be changed during the operational lifetime at the 
end of a groundtrack repeat interval. 

Determination of the satellites’ orbits will be provided by one of three methods, 
all of which will have almost complete coverage of the orbits. The orbit tracking will 
use the Tracking Data Relay Satellite System (TDRSS), the Global Positioning System 
(GPS), and/or the ground based doppler tracking network, TRANET. The system that 
is finally selected must be able to support the required orbit accuracy. The 3a accuracy 
for the gravitational part of the mission has been specified as 100 m in the radial 
direction and 300 m in the along- and cross-track directions. The three tracking 
systems under consideration have precisions that readily support these orbit 
requirements. For the magnetic mission, the orbits must be known to within 60 m for 
radial, and 100 m for the along- and cross-track directions [Keating, et al., 7986]. 
The criteria for the gravitational aspect of the mission are discussed in more detail in 
Chapter 4. 
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2.5 Specifications for Mission Simulation 

Specific mission characteristics have not been chosen therefore, nominal 
mission characteristics have been selected for this study that reasonably represent the 
actual mission. The gravity parameter, (I, used in the simulation for this study equals 

r 'i o 

3.9860064 x 10 km /sec and the mean equatorial radius was taken to be to 6378. 145 
km. Although the semimajor axis will vary slightly depending on the geopotential field 
used, it will maintain a mean height above the reference ellipsoid of approximately 160 
kilometers. 

The choice of a 160 km altitude produces exactly 525 revolutions of the 
satellites in 32 sidereal days. Since 525 and 32 are not commensurate, the orbits are 
guaranteed to have their first repeat in exactly 32 sidereal days [Thobe and Bose, 
1985]. A sufficiently long time interval for the repeat of the groundtrack must be used 
in order to resolve the order of the spherical harmonic expansion of the geopotential 
field. The number of revolutions in a groundtrack repeat interval is twice the highest 
harmonic order that can be determined, and thus the recovery of the geopotential 
coefficients up to degree and order 262 is theoretically possible [Colombo, 1984]. 

The model used in this study, except for Chapter 6, included the Earth’s gravity 
field as the only force acting on the satellites. No luni-solar, external gravitational 
forces or nongravitational forces were considered. The coordinate system was body- 
fixed with constant Earth rotation, i.e. precession, nutation, and polar motion were not 
considered in the model. 
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CHAPTER 3 


FROZEN ORBITS 


3.1 Background on Frozen Orbits 

For satellite science experiments which require the same point on the Earth to be 
sampled numerous times, a repeating groundtrack is necessary. In addition, if it is 
desired that along a particular latitude the altitude variations have a constant mean value, 
then a "frozen orbit" is required for this application [McClain, 1987]. Both of these 
characteristics are employed in the recovery of the gravity field harmonics for the 
Geopotential Research Mission. A frozen orbit's shape is held constant, thereby 
minimizing the variations in the altitude; in addition, a frozen orbit helps to maintain the 
nature of the groundtrack repeat [Nickerson, Herder, Glass, and Cooley, 1978], 

Several satellite missions have used the frozen orbit concept, including the 
Atmospheric Experiment satellites (AE-3 and AE-5), SEASAT, LANDSAT, the Heat 
Capacity Mapping Mission (HCMM) [Nickerson, et al., 1978] and GEOS AT. Some 
of the results from studies of these missions are discussed in this and in later chapters. 

For the Geopotential Research Mission, as well as some other missions, an 
exact groundtrack repeat is required. Rotation of the line of apsides due to even zonal 
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harmonics makes an exact repeat more complicated to achieve. At critical inclination, 
perigee does not precess; however, because of mission considerations and constraints, 
not all missions are allowed to have this particular inclination of 63.4°. To facilitate the 
accomplishment of an exact repeat, it would be convenient if the perigee motion could 
be prevented at noncritical inclinations. For low Earth satellites, the line of apsides will 
undergo a secular motion of approximately four degrees per day, considering only the 
J 2 effect [Roy, 1978]. Therefore, to obtain a frozen orbit, other gravity perturbations 
must interact with the secular effects caused by the even zonal harmonics to produce 
small or zero motion in the argument of perigee for noncritical inclinations. Using a 
frozen orbit, where perigee will not precess, enables an exact groundtrack repeat to be 
obtained more easily. 

If Fourier transforms are used as the technique to recover the gravity field, then 
constant mean altitude over the subsatellite points is needed [Colombo, 1985]. Since 
for the GRM satellite orbits only gravitational perturbations are considered, then the 
semimajor axis, the eccentricity and the inclination will have constant mean values. 
However, the long period motion in eccentricity and the secular rate in the argument of 
perigee need to be removed in order to properly employ the Fourier transform 
technique. This can be achieved with a frozen orbit A frozen orbit can also maintain a 
lower mean eccentricity than a non-frozen orbit [Nickerson, et al., 1978]. This chapter 
investigates the need for frozen orbits, demonstrates the manner in which frozen orbits 
are derived, and describes the characteristics of frozen orbits. 
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3.2 Definition of Frozen Orbit 

The definition of a frozen orbit, according to Kozai [1959], is an orbit that has 
no secular precession of the argument of perigee at any inclination; the longitude of the 
ascending node is not required to be constant. In addition to eliminating the precession, 
the argument of perigee and the eccentricity of the frozen orbit will not display any long 
period effects due to the odd degree zonal harmonics. 

Cook [1966] determined that if only the long period and secular effects of a 
disturbing function are considered, then a value for eccentricity can be found for a 
given semimajor axis and inclination that will eliminate the unbounded motion of the 
line of apsides for nearly circular orbits. The orbit will remain frozen as long as 
nongravitational perturbations do not interfere with the orbital motion. Cook 
demonstrated this result with a disturbing function that included only one even degree 
harmonic (J 2 ), and all the odd degree zonal harmonics to J 9 . 

Since the argument of perigee is not well defined for near-circular orbits, the 
rate of change of the argument of perigee will not have a smooth secular change, but 
instead, it will exhibit nonlinear variations. Such nonlinear effects allow the long 
period trends due to the odd zonal harmonics to effectively cancel the secular trends due 
to the even zonal harmonics. In order to obtain a frozen orbit, the geopotential field 
used to calculate the orbital elements must have the minimum of one odd and one even 
degree zonal coefficients, e.g., J 2 and J 3 . 


The addition of J 4 to Cook's analysis changes the value for the frozen orbit 
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eccentricity by approximately 3%. As more even zonal harmonics of higher degree are 
added for the calculation of the frozen orbit eccentricity, the changes in the eccentricity 
become less significant [Nickerson, et al ., 1978] because of the general reduction in 
coefficient magnitude with the increase in degree. The value of the eccentricity, 
calculated for a geopotential field of even and odd zonal harmonics to J 9 , will require 
little adjustment if a much larger geopotential field is used, including the tesseral 
harmonics. Although Cook did not include any of the even zonal harmonics above J 2 
in his work, the calculations for the frozen orbital elements required for this study 
included J 4 . 

Lagrange's planetary equations for the Keplerian orbital elements contain a 
singularity when the eccentricity equals zero. To remove the singularity, and to have a 
set of orbital elements that will be well behaved for small eccentricities, a modified set 
of orbital elements was chosen. The elements £, r|, and a replaced e, CD, and t p in the 
standard, Keplerian set of orbital elements. 

£ = e cos co 
T| = e sin co 
ct = co + nt p 

where n is the two-body mean motion, t p is time of perigee passage, e is the 
eccentricity, and co is the argument of perigee. Semimajor axis, inclination and 
longitude of the ascending node complete the set of elements. Lagrange's planetary 
equations for the modified orbital elements are provided in the Appendix A [Taft, 
1978]. The two orbital elements, T] and £, are the components of the eccentricity 
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vector, which points in the direction of periapsis from the center of the coordinate 
system [Bate, Mueller and White, 1971]. 

By neglecting the short period effects, an analytical solution for the time 
derivatives of T| and £, can be used instead of actually integrating the equations of 
motion provided in Appendix A. The solutions are: 


% = A cos ( kt + a ) 
rj = A sin ( kt + a ) + C/k 


(3.1) 


The equations for C and k are also provided in the Appendix A, as given by Cook 
[1966]. The secular effects due to the even zonal harmonics are contained in k and the 
long period effects due to the odd zonal harmonics are contained in C. The amplitude, 
A, for the oscillation in % and T| is independent of the amplitude of the long period 
effects produced by the disturbing function. The amplitude and phase angle, a, depend 
only upon the initial conditions. 

3.3 Characteristics of Frozen Orbits 

When the orbit motion is plotted in the ( "q ) plane, two possibilities can 
occur. First, if A > I C/k I, then the argument of perigee will be unbounded, and the 
orbit will not be frozen (Figure 3.1a). The second possibility is that if A < I C/kl, in 
which the argument of perigee will be bounded between 0° and 180°, and will oscillate 
around 90° as illustrated in Figure 3. lb [Cook, 1966]. When the eccentricity is exactly 
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equal to C/k, then the amplitude will be zero and the argument of perigee will be held 
fixed at 90° instead of oscillating, and the orbit is termed "frozen" [Cook, 1966]. 

The time required to complete the circles in Figures 3.1 is equivalent to the 
longest period in the disturbing function due to the odd degree zonal harmonics; in the 
case of GRM, this interval is about 79 days. It should be noted that ri oscillates about 
C/k in Figure 3.1b, but maintains a value greater than zero; when this is the case, the 
orbit is termed a "frozen" orbit. The eccentricity will be constant when it is equal to 
C/k, and it will not exhibit long period perturbations due to the odd zonal harmonics as 
a function of time. If the eccentricity is not exactly equal to C/k, then a long period 
oscillation will occur, but the amplitude of that oscillation depends only upon the initial 
j conditions of Equations (3.1) [Cook, 1966]. 

[ The mission simulation selected for this study of GRM has 525 exact 

revolutions in 32 sidereal days, which requires that the semimajor axis equals 6526.988 
| km, and the eccentricity equals 0.00153084 for the disturbing function used by Cook. 

For the inclination equal to 90° with conditions given above, the phase plane diagram 
and the eccentricity as a function of the argument of perigee are displayed in Figures 
3.2a and 3.2b. The center represents the frozen orbit (e 0 = 0.00153084) and was 
based on the even degree zonal harmonics, J 2 and J 4 , and on all the odd degree zonal 
harmonics to J 9 . The amplitude. A, is equal to C/k - e Q and the phase angle, oc, is set to 
90°. In contrast, the non-ffozen orbit figures for the case where A > I C/k I are 
provided in Figures 3.3a and 3.3b. 
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From Kaula [1966], the deviations from a secularly precessing orbit due to the 
disturbing function for eccentricity and argument of perigee are: 

= iiae 1 F 1mD G| na (l-g 2 ) 1/2 [(l-g 2 ) 1/2 (l-2p+q)-(l-2p)]Simpq 
" na 1+ V[(l-2p)d> + (l-2p+q)A* + m(Q-©)] 

AcOtapq = Mae 1 [(l-g 2 ) 1/2 Fi mp QG| pq /ag ) - cot/ (l-e 2 y m (dF lmr Jdi )Gi pq ]S~ lmpq 
n a + e [(l-2p)d) + (l-2p+q)Aif + m(Q-0)] 

where ae is the mean equatorial radius of the Earth, © is the rotation rate of the Earth, 
Gipq is the eccentricity function, Fi mp is the inclination function, and Si mp q is the 
integral of S]jn P q with respect to its argument. Since only zonal harmonics are being 
considered: 

Simpq = Qm cos[(l-2p)co + (l-2p+q)Af ] + Cim sin[(l-2p)co + (l-2p+q)M ] 

When the denominators of Equations (3.2) become exactly zero, as in the case 
of zero eccentricity or in the case of resonance (Simpq -0), these equations become 
invalid [Kaula, 1966]. The long period effects generated by odd zonal harmonics are 
associated with the coefficient of the time derivative of the argument of periapsis, cb. 
These terms will be zero for the frozen orbit, and since only zonal harmonics are 
considered, m is zero. The remaining term is the time derivative of the mean anomaly 
term, M, which is associated with the short period terms. Cook's theory only 
considers the secular and the long period effects and neglects the effect of the short 
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period terms, therefore, the entire denominators of both Equations (3.2) are zero, and 
Kaula's equations are not valid to use for Cook's frozen orbits. 

3.4 Addition of Short Period Terms 

With the inclusion of the short period perturbations, the closed form of the 
solutions for | and f| can no longer be used. Instead of the zonal harmonics to J9, the 
full 9x9 geopotential field was used to create a frozen orbit that, unlike the preceding 
analysis, includes short period effects. The zonal harmonics and the 9 x 9 geopotential 
field were taken from the OSU322 field described in Section 1.2. The orbit was 
generated using the software UTOPIA, the University of Texas Orbit Processor; which 
will be described briefly here, but in more detail in Chapters 5 and 6. The numerical 
integration of the equations of motion with UTOPIA was carried out for 32 sidereal 
days. 


The frozen orbit conditions were met by choosing the mean orbital elements to 
be the frozen orbital element values, arrived at through Cook’s equations, as inputs to 
SPENEW. The software SPENEW evaluates the analytical expressions for a secularly 
precessing ellipse and generates a position ephemeris file for 32 days. Points from the 
secularly precessing ellipse for the frozen orbit were used as the observations for 
UTOPIA. A least squares fit to the position ephemeris was made by UTOPIA to obtain 
a set of initial conditions that should remain frozen in a mean sense, and to study the 
short period effects on the frozen orbit due to the gravity harmonics. 
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The specific characteristics of the file created by SPENEW were as follows: 

a = 6526988 m 
e = 0.00153084 
i = 90° 

co = 90° d) = 0 rad/sec 

Q = 90° Q = 0 rad/sec 

M = 0° M- 0.1 19636321 x lO 2 rad/sec 

The results from the least squares fit by UTOPIA to the ephemeris generated by 
SPENEW are shown in Figures 3.4 and 3.5 which now include the short period 
effects. From the phase diagram for the non-frozen orbit (Figure 3.1a), the value for T| 
becomes less than zero, therefore, the argument of perigee will not be bounded and the 
orbit cannot be frozen. This appears to occur when the short period effects are included 
(Figure 3.4a and b). One complete trace of the curves in these figures is performed in 
one orbital period. The patterns are due to the inclusion of the short period effects only 
and are not associated with the long period tracings seen in Figures 3.2a and 3.2b. The 
amplitude of the short period effects are sufficiently large to apparently destroy the 
frozen orbit integrity. However, the changes in the orbital elements with respect to time 
(Figures 3.5a and 3.5b) indicate that their mean values remain frozen. Consequently, 
even though the short period terms increase the osculating values of the eccentricity and 
the argument of perigee beyond the bounds of the frozen orbit, the mean orbit elements 
retain the frozen characteristics. 


i 
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From the individual orbit element plots (Figure 3.5a and b), the mean values of 
argument of perigee and eccentricity remain constant at approximately their frozen orbit 
values. There is no long period trend in either of the elements, and there is no secular 
trend in argument of perigee. Though the osculating to varies from 0° to 360°, the 
mean value remains constant at 90°» allowing a frozen orbit to still exist even in the 
presence of short period effects [Nickerson, et al., 1978]. This effect will be illustrated 
more clearly in later sections when the trends in the orbit element plots for a frozen orbit 
are compared to the orbit element characteristics of the non-frozen orbit. 

3.5 Non-frozen Orbits 

As a comparison to the frozen orbit, a non-frozen orbit was generated using the 
same 9x9 geopotential field used to generate the frozen orbit, described in the 
preceding section. However, for the non-frozen orbit, d> was equal to -4.59 
degrees/day, instead of the frozen orbit value of zero. The time rate of change of the 
argument of perigee, as well as the time rate of change of the mean anomaly, must be 
included when determining the value for the semimajor axis that will provide an exact 
groundtrack repeat for the non-frozen orbit case. 

If the orbit is frozen, d) = 0, then only M must be considered in the calculation. 
As discussed earlier, the satellites for GRM in this study must repeat their groundtracks 
every 32 sidereal days. With the inclusion of a nonzero d>, the periods of the orbits will 
not be changed, but the value for the semimajor axis for a non-frozen orbit will differ 
from the frozen orbit's value. The equation for d), provided by Kaula [1966], 
calculated as a function of J 2 and J 4 only is as follows: 
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d) = -3nJ 2 a P 2 (5sin 2 / -4) - 5nJ 4 a P 4 ( 105/64 sin 4 / - 15/8 sin 2 / + 3/8) 

vi/i 2x2 2 /-% 2\4 4 

4(l-e ) a (l-e ) a 

The time rate of change of mean anomaly plus the time rate of change of 
argument of perigee must sum to the value of the mean motion that will contain 525 
revolutions in 32 sidereal days. The resulting mean semimajor axis for the non-frozen 
orbit is 6523.608 km. The angular rates are: d) = -0.928769168868 x 10 -6 rad/sec. Cl 
is still equal to zero, and M = 0.119729152669167 x 10 -2 rad/sec. Figures 3.6a and 
3.6b are the orbit element plots for the non-frozen orbit for 64 days, two complete 
groundtrack cycles. A long period trend is apparent in the eccentricity plot, which has a 
value of about 79 days. The argument of perigee has an obviously secular trend, unlike 
the frozen orbit co. Inclination and longitude of the ascending node plots are not given 
since both of these parameters remain essentially constant for the polar orbit. 

The initial latitude and longitude for the leading GRM satellite used in this study 
were 88.688° and 169.757°. After integrating for 32 sidereal days, with a 9 x 9 
geopotential field, the final longitude and latitude for the non-frozen orbit were 88.689° 
and 169.757°. After 64 sidereal days, the latitude and longitude for the satellite were 
88.562° and 169.749°. After 32 sidereal days, the non-frozen orbit satellite had a 
nearly exact groundtrack repeat but, after 64 sidereal days, the deviation in latitude was 
0.126°. For the frozen orbit with a 9 x 9 geopotential field, the deviation in latitude 
after 64 days was only 0.005° (which will be demonstrated in Section 5.2), indicating 
the greater ease in maintaining a repeating groundtrack for a frozen orbit than for a non- 
frozen orbit. Further comparisons of the frozen and non-frozen features, in particular 
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the drift problem between the two satellites, are presented in Chapter 4. 

The phase plane and eccentricity versus the argument of perigee diagram for the 
non-frozen orbit are illustrated in Figures 3.7a and 3.7b. The satellite's orbit was 
integrated for 80 days in order to complete one cycle of the longest period due to the 
odd zonal harmonics. The frozen orbit, which includes short period effects, has a very 
distinct pattern, with clear borders (Figures 3.4a and 3.4b). The non-frozen orbit, 
however, produces a more diffuse pattern. Also, the magnitude of the parameters are 
larger than the frozen orbit's values, indicating larger variations in the orbit element 
values. 

3.6 Summary 

This aspect of the study was performed using a 9 x 9 geopotential field. A 
study presented by Schutz, Tapley, Lundberg and Halamek [1986] contains the phase 
plane diagrams for a 180 x 180 geopotential field. The results from the 180 x 180 field 
are nearly identical to the phase plane diagrams presented here (Figures 3.4a and 3.4b), 
thereby suggesting that the dominant short period effects are contained within the 9 x 9 
field, and that the dominant secular and long period effects are due to the first nine 
zonal harmonics. The short period amplitudes due to h alone are at least three orders 
of magnitude larger than the short period amplitudes due to any of the other harmonic 
terms [Kaula, 1966]. 

As was explained in this chapter, the importance of frozen orbits for the 
Geopotential Research Mission is twofold; frozen orbits allow repeating groundtracks 
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to be more easily maintained and they provide a near constant altitude over individual 
subsatellite points. Since it is desired to repeat the entire mission several times to 
secure reliable measurement data, the repeating groundtrack is a necessary mission 
requirement. For certain geopotential recovery techniques minimizing altitude 
variations is essential, therefore, a frozen orbit will be required to meet this condition as 
well. 


Results in Chapter 4 will demonstrate that a frozen orbit allows the repeating 
groundtrack to be more easily maintained than a non-frozen orbit over the mission 
lifetime. Results in Chapter 6 will provide an indication of the frozen orbit stability in 
the presence of other gravitational perturbations. 
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PHASE PLANE PLOTS FOR FROZEN ORBITS 
NO SHORT PERIODIC TERMS INCLUDED 
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Figure 3.2a 

Phase plane for frozen orbits 





ECCENTRICITY *10 


40 


VARIATIONS OF ECCENTRICITY WITH ARGUMENT 
OF PERIGEE FROZEN ORBITS 


N 

0 > 



ARGUMENT OF PERIGEE (DEGREES) 


Figure 3.2b 

Variation of eccentricity versus argument of perigee. 
No short period effects were included: Frozen orbit 
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Figure 3.4a 

Phase plane plot: Short period effects were included. 
Mean orbit elements were frozen. 
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ECCENTRICITY VERSUS ARGUMENT OF PERIGEE 
FROZEN ORBIT 9X9 GEOPOTENTIAL 
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Figure 3.4b 

Eccentricity versus argument of perigee: Short period effects were included. 
Mean orbit elements were frozen. 
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FROZEN ORBIT IN A 9 X 9 GEOPOTENTIAL 
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Figure 3.6a 

Eccentricity versus time 
For 64 sidereal days: Non-frozen orbit 
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PHASE PLANE FOR NON-FROZEN ORBIT 
GEOPOTENTIAL FIELD IS 9 X 9 
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Figure 3.7a 

Phase plane for non-frozen orbit 
Short period effects were included 
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ECCENTRICITY VERSUS ARGUMENT OF PERIGEE 
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Figure 3.7b 

Variation of eccentricity versus argument of perigee 
Short period effects were included 
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CHAPTER 4 


ADJUSTMENT OF INITIAL CONDITIONS 


4.1 Introduction 

A set of initial conditions must be determined for both satellites that will best 
satisfy the mission requirements. Since it is expected that the trailing satellite follows 
the lead satellite in nearly the same orbit, only one set of orbital elements for one of the 
satellites will be discussed. 

As described previously in Section 2.5, the satellites in this simulation will have 
a 32 sidereal day repeat with 525 exact revolutions in that time. The orbital period that 
will produce these conditions can be determined by dividing the number of seconds in 
32 sidereal days by 525. This orbital period produces a resulting mean motion of 
1.1963627 -x 10' rad/sec. If the orbit is not polar, the determination of the period is 
calculated by: 


P = (l+Q/co e ) ND/NR 

where P is the orbital period, £2 is the mean time rate of change of longitude of the 
ascending nodes, co e is the rotation rate of the Earth, ND is the integer number of days 


of the groundtrack repeat interval, and NR is the integer number of revolutions in the 
repeat cycle. If the orbit is polar, as the case of GRM, the precessional rate of the node 
will equal zero and the orbital period will simply be equal to the ratio ND/NR. 

The time rate of change of the mean anomaly is equated to the value of the 
mean motion given above. However, the secular time rate of change of mean anomaly 
consists of not only the two-body term, but it is also a function of the even zonal 
harmonics. The equation for the time rate of change of the mean anomaly, til, as a 
function of the disturbing terms due J 2 and J 4 only is: 

til = n [ 1 + 3J 2 a e 2 /{ 4(1 -e 2 f 2 a 2 } (3cos 2 i - 1) + 3J 4 (a Ja ) 4 (105/64 sin H - 
15/8 sin 2 / + 3/8) ( 1-e 2 )' 5/2 { 1 - (1 + 3/2e 2 )/(l - e 2 )} ] 

where n is the two-body mean motion. It should be noted that the time rate of change 
of mean anomaly is a function of semimajor axis and eccentricity. In order for the 
satellites to repeat in 32 sidereal days, the value of the semimajor axis, a, must satisfy 
the time rate of change of mean anomaly when equated to the previously determined 
value for the mean motion of the orbit (calculated from the ratio 2% NR/ND). The 
eccentricity must satisfy the frozen orbit conditions defined earlier. 

A frozen orbit is created by selecting the mean orbit elements to be a = 
6526988.0 m, e = 0.001534965, i = 90°, co = 90°, Q = 90°, and M = 0°. The time 
rate of change of argument of perigee is zero, complying with the conditions for a 
frozen orbit. A least squares fit to this frozen orbit trajectory was made to calculate the 
initial conditions that best fit the frozen orbit trajectory for 32 sidereal days. The 
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generation of this mean, frozen trajectory was described in Section 3.4. The initial 
conditions will depend upon the size of the geopotential field used, the gravity 
parameter, |i, and the mean equatorial radius, ae. The gravity fields used in this chapter 
are taken from the OSU322 field [Rapp, 1981] and from OSU86F [Rapp and Cruz, 
1986] described in Section 1.2. 

Colombo [1985] provided initial conditions for a geopotential field that 
consisted of the OSU322 36 x 36 subfield plus the zonal harmonics to degree 300, with 
no other temporal gravitational effects (e.g. tides, precession, etc.). These initial 
conditions used a gravity parameter of 3.986013 x 10 km /sec and a mean equatorial 
radius of 6378.155 km. Colombo's initial conditions in Cartesian coordinates are 
provided in Table 4.1. 


Table 4.1 


Colombo's Initial Conditions for Leading and Trailing Satellites 




Position (ml 


Velocity (fn/s'l 

Leading Satellite 

r x 

0.0 

v x 

0.0 


r y 

-150000. 

Vy 

-7817.1468749596 


r z 

6514763.771448 

v z 

-179.513061923 

Trailing Satellite 

r x 

0.0 

V X 

0.0 


r y 

150000. 

v y 

-7817.1468749596 


r z 

6514766.990466 

v z 

179.513061923 
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5 3 2 

Using a more current gravity parameter, 1 1, of 3.9860064 x 10 km /sec and a 
mean equatorial radius of 6378.145 km, a new set of initial conditions were determined 
for a 36 x 36 OSU322 subfield which will be described in Section 4.8. 

4.2 Determination of a Principal Set of Initial Conditions 

A 32 sidereal day repeating groundtrack trajectory was obtained by generating a 
secularly precessing elliptical orbit defined by the frozen orbit conditions provided in 
the previous section. As described in Section 3.4, UTOPIA was used to make a least 
squares fit to this frozen orbit. The initial osculating position and velocity for a 
geopotential field with only zonal harmonics to J 9 are presented in Table 4.2. 


Table 4.2 

Initial Conditions for Zonal Harmonics to J 9 

Position (m) 


Velocity (m/s) 

0.0 

v x 

0.0 

-32010.44187981 

Vy 

-7819.018360577 

6516497.39025295 

V Z 

-37.86736549955 


These initial conditions were numerically integrated forward and backward in time until 
y equaled ±150000 meters, producing initial conditions for both the leading satellite and 
the trailing satellite. These y-coordinates were selected since they allow the initial 
conditions to resemble Colombo's initial conditions, and caused the satellites to begin 



at approximately 300 km apart; the actual separation distance will be a mission 
specification. The epoch time used by Colombo and assumed for this study was the 
Julian date of 2445700.5, January 1, 1984, 00:00. 

The initial conditions for both satellites, with a geopotential consisting of zonal 
harmonics to J 9 only, were numerically integrated for 32 sidereal days to determine 
their relative range, relative range-rate, and relative acceleration between the satellites. 
The relative range values were subtracted from a 300 km separation distance to assess 
the ability of the initial conditions to produce little or no secular trend in the relative 
motion. If a secular trend (or drift) does appear, it can be eliminated by adjusting the 
initial conditions. A positive secular trend indicates that the satellites are drifting 
together, and they are drifting apart if the trend is negative. Both of the satellites are 
required to have an exact repeat after 32 sidereal days, that is, the longitude and latitude 
of both satellites must return to their initial values. If only zonal harmonics are 
considered, there will be no longitudinal error, since the two orbits will remain in the y- 
z plane. However, when the tesseral and sectorial harmonics are included in the 
geopotential model, longitudinal errors in the groundtrack will be introduced. 

If the initial conditions created for the case with zonal harmonics only are used 
in a numerical integration with a 9 x 9 geopotential field, there will be a resulting drift 
generated between the two satellites, and their final subsatellite points after 32 sidereal 
days will not be equal to their initial values; therefore, they will not produce a repeating 
groundtrack. Instead of fitting the ephemeris created from a 9 x 9 geopotential (or any 
other geopotential for that matter) to the frozen orbit defined by a secularly precessing 
ellipse, an adjustment can be made in the initial conditions derived for the zonal 
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harmonics only geopotential field to each satellite, based on the fact that the two 
satellites must repeat after 32 sidereal days. 

4.3 Using 32 Sidereal Days to Determine Initial Conditions 

A linear approximation for the required change in the radial direction is used to 
simultaneously eliminate the drift rate and the latitude errors. Since the orbits are polar 
and the initial conditions are close to the north pole and the orbits are nearly circular, 
then a, r, and z are approximately equal. As a simplification, the corrections could be 
applied in the z direction, however, if the satellites are located somewhere other than 
above the pole, the correction should be made in the radial direction. The change in 
mean motion, where two-body mean motion only is considered, is equal to twice the 
drift rate between the two satellites. Taking the partial derivative of the two-body mean 
motion, n, with respect to the mean semimajor axis, a: 

1/2 DR = 5n = -3/2 (\yja s ) m Azo (4.1) 

where DR is the drift rate between the satellites and Azq is the total change needed in the 
initial conditions to correct for the drift between the satellites. From Equation (4.1), the 
net change in the z-coordinate for the satellite position that will eliminate the drift rate 
can be determined. Note, for consistent units, the radius, r, must be included in the 
right-hand side of the Equation (4.1). 

With Azq due to drift, the corresponding latitude change that will be incurred 
can be calculated. Since the orbit is near circular, the ratio of the latitude change due to 
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drift (A<{>, which is equal to 5 rt) to the change in initial z to compensate for the drift only 
(Azq) remains essentially constant. After a 32 sidereal day integration of both satellites, 
the latitude difference (5 <{>,•) between the initial conditions and the final conditions of 
each of the satellites can be computed. The error due to the drift is corrected by 
assigning Azq to the satellite that has the largest individual latitude error. The excess 
latitude difference is then adjusted by using the ratio of the change in latitude to the 
change in z, A^/Azq, due to drift: 

( A(J) / Az 0 ) Drift = ( 8<J>; /8z 0 ; ) (4.2) 

where 8 zq ,• is the required change in the initial z component to obtain an exact repeat in 
latitude for that particular satellite, and i indicates the satellite being considered. The 
program FIXDRF, given in Appendix B, was used to calculate the required corrections 
in each satellites’ initial conditions. 

Note that the 8z 0 i correction will be the same for both satellites; this prevents 
inadvertently generating drift between the satellites at the possible sacrifice of the orbit 
of one satellite not repeating as well as the other. However, the total amount of the 
correction in z added to each satellites' positions will be different, since the satellite 
with the largest change in latitude will also include the adjustment for the drift rate. The 
determination of whether total correction for the z component is added (or subtracted) 
from the initial conditions, is dependent on whether the satellite has moved ahead (or 
has trailed behind) at the final time from its latitude at the initial time. If the satellite is 
behind, then the correction must be subtracted from the initial conditions to increase the 
speed of the satellite (Figure 4.1a). If it is ahead, then the correction must be added to 


58 


raise the satellite's altitude and slow it down (Figure 4.1b). 

Small changes in z affect latitude, but have little effect on longitude. To correct 
the error in longitude, which will occur when sectorial and tesseral harmonics are 
included in the geopotential, the x-coordinate must be adjusted. Small changes in x will 
affect the longitude, but not significantly change the latitude (Figure 4.2a). Changes in 
the y-coordinate have a similar effect as changing the z-coordinate, that is, a deviation 
in y adjusts latitude (Figure 4.2b). Without precession and nutation, the offset in 
longitude is the same for both satellites and, therefore, the correction to the x 
component of the initial conditions will be the same for both satellites. 

For longitude, unlike in the latitude adjustment, when the initial value for x is 
changed, the initial longitude value also changes. This makes an analytical approach 
for finding the needed x adjustment (8x 0 ) difficult to obtain. Several different 
adjustments in the initial x-coordinate were made using the 9 x 9 geopotential field. 
The other components were held fixed, as the latitude error has already been corrected. 
The resulting changes in longitude are plotted against the given changes in x, and a least 
squares quadratic fit was made through the points to determine the 8x 0 that corresponds 
to 8 A, = 0° (Figure 4.3). It was found that the same value for 8x 0 can be used to 
correct the longitude offset in the satellites' initial conditions for a 36 x 36 geopotential 
field as was needed for the 9x9 geopotential field. This implies that the dominant 
terms that cause the longitude offset are contained within the 9 x 9 geopotential field 
and that the higher degree and order terms have little affect on longitude changes (i.e., 
the orbital period does not change significantly). 
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The initial conditions in a 9 x 9 geopotential field taken from the OSU322 field 
that will have a 32 sidereal day repeat both in geodetic latitude and in longitude (to three 
decimal places) are presented in Table 4.3a. The resulting latitudes and longitudes for 
each of the satellites are also presented (Table 4.3b). These results indicate that the 
satellites will have a closure of within 450 meters, even after 64 sidereal days (two 
complete repeat groundtrack cycles). 


Table 4.3a 

Initial and Final Conditions for the 9 x 9 Simulation 
Repeating groundtrack 




Position (ml 


Velocitv(m/s1 

Leading Satellite 

rx 

253.7524 

v x 

0.0 


r y 

-150000. 

v y 

-7816.570937552 


r z 

6515233.9274476 

v z 

-179.497721338 

Trailing Satellite 

r x 

253.7524 

v x 

0.0 


r y 

150000. 

Vy 

-7816.570937516 


r z 

6515238.3234627 

V Z 

179.5005234656 
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Table 4.3b 
Frozen Orbits 

Latitude and Longitude Values for each Satellite 

Initial Lead Satellite Trailing Satellite 

Latitude 88.690° 88.690° 

Longitude 169.757° 349.563° 

Final after 32 Sidereal Davs 

Latitude 88.691° 88.690° 

Longitude 169.757° 349.564° 

Final after 64 Sidereal Davs 

Latitude 88.684° 88.694° 

Longitude 169.757° 349.563° 


The drift rate between the two satellites in the 9 x 9 geopotential field with the given 
initial conditions was 0.15564 m/day. Figures 4.4a, 4.4b, and 4.4c illustrate the 
relative range, range-rate and acceleration for these initial conditions. Note that the 
secular trend between the frozen orbit pair is small, indicating little drift between the 
satellites (Figure 4.4a). 

The relative motion figures for the non-frozen case, described in Section 3.4 
with the 9 x 9 geopotential field, are presented in this section as a comparison to the 
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frozen orbit relative motion. Even though the satellites were close to an exact repeat for 
the first 32 days, a drift between the two non-ffozen orbiting satellites of approximately 
-50 m/day results. If the z component of position for the leading satellite is adjusted 
slightly, the drift can be eliminated, but the satellite will no longer repeat as accurately. 
Figures 4.5a, 4.5b and 4.5c present the relative motion after 64 sidereal days, and the 
latitudes and longitudes of each satellite are provided in Table 4.3c. Clearly, a less 
accurate repeating groundtrack results for the two, non-frozen satellites compared to the 
frozen orbit. 


Table 4.3c 
Non-ffozen Orbits 

Latitude and Longitude Values for each Satellite 


Initial 

Lead Satellite 

Trailing Satellite 

Latitude 

88.688° 

88.688° 

Longitude 

169.757° 

349.563° 

Final after 32 Sidereal Davs 


Latitude 

88.689° 

88.688° 

Longitude 

169.756° 

349.563° 

Final after 64 Sidereal Davs 


Latitude 

88.814° 

88.797° 

Longitude 

169.785° 

349.533° 
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4.4 Non-Polar Adjustments 

For this study, the satellites were initially located near the poles, therefore, the 
change in radius was approximately equal to the change required in the z-coordinate. If 
the initial conditions are not near the pole, then the corrections would have be added 
radially to the orbit. Since the orbit plane is mainly in the y-z plane, the x-coordinate 
can be neglected. 

(5r c ) 2 = (5y 0 ) 2 + (Sz 0 ) 2 

5zo / 8y 0 = tan(<J) 0 ) 

where <j> 0 is the initial latitude of the design orbit. The required change in radius is 
determined by the previously described method. With two equations and two 
unknowns, the change in initial y- and z-coordinates can be calculated regardless of the 
initial latitude. 

4.5 Sensitivity Study 

With the initial conditions from the 9 x 9 geopotential termed "nominal", a 
sensitivity study was conducted. Table 4.4 contains the results of several cases that 
investigated the sensitivity of the satellites to perturbations in initial position ranging 
from one centimeter to ten meters. These perturbations were made to determine the 
sensitivity of the relative motion to errors in the initial conditions in order to observe the 
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corresponding changes in the groundtrack of the satellites. It should be noted that a 
change in z, at least for a change no larger than 10 meters, will not change the initial 
values for the latitude and longitude of the altered satellite. Once a satellite is ahead, or 
behind, of its designed trajectory, it will remain ahead, or behind. Therefore, the cases 
were divided into three groups. For the first two groups, the leading satellite remained 
behind its designed orbit, because the perturbation is in the positive z direction. For the 
third group, the leading satellite was ahead of its designed orbit. 

The first group had a perturbation in the positive z direction for the lead satellite 
only, and the second, trailing satellite remains unperturbed. Four deviations from the 
nominal were made: +1 cm, +10 cm, +1 m and +10 m, in the positive z direction for 
the leading satellite only. 

For the +1 cm case, the final conditions for latitude and longitude were equal to 
the final conditions for the nominal case, to at least three decimal places and the 
groundtrack maintained the of closure of within 100 m. A drift rate between the two 
satellites of 3.24926 m/day was acquired, compared to 0.155642 m/day from the 
nominal case. This drift rate is very small and would bring the two satellites closer 
together by only 100 meters in 32 days. 

For the +10 cm case, the final conditions for latitude and longitude errors were 
0.009° and 0.001°. The groundtrack difference after 32 days has increased an order of 
magnitude over the +1 cm case to 1000 m. An acceptable error in the groundtrack 
closure is considered to be 10 km [Schutz, et al., 1986]. The drift rate also increased 
an order of magnitude to 31.123122 m/day. This would bring the two satellites about 
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1000 meters closer than their initial 300 km separation distance. 

For the +1 m case, the errors in latitude and longitude were 0.082° and 0.007°. 
The drift rate increased another order of magnitude to 309.852756 m/day. The 
difference in the groundtracks between this and the nominal case after 32 days is just 
under 10 kilometers. It should be noted that this case brings the groundtrack errors 
close to the limitations that there must be a 10 km closure for the groundtrack repeat. A 
one meter error in the initial conditions is the upper limit allowable, based on the 
analysis presented here. 

The last case for this group was a +10 m perturbation which resulted in latitude 
and longitude errors of 0.864° and 0.189°. The drift rate was 3097.8527 m/day. This 
would result in the satellites being roughly 200 km apart after thirty-two days which is 
one hundred kilometers less than the mission specification. The groundtrack had a 
closure error of 100 kilometers when compared to the nominal. From these results, an 
error of only 10 meters in one of the satellites would produce unacceptable orbit 
conditions and would not meet the mission specifications. 

The second group had perturbations added to both satellites’ initial conditions in 
the same direction; positive z. Both satellites were traveling behind their designed 
orbits, but since the perturbations were equivalent for both satellites, they produced no 
net drift rate between them because they are traveling at the same relative rate. The 
errors in longitude and latitude remained the same as they were for the leading satellite 
from the first group (Table 4.4). 



T iM e 4.4 

Perturbation in Initial Conditions tor l ,aulin<? and Trailing Satellites 
and the Resulting Errors in the Grouin!i;: , Hc Repeat 
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-10 m / +10 m -0.864 0.038 96.2 0.865 0.039 96.2 6183. 
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The third group had equivalent perturbations, but they were applied to the 
satellites in the opposite direction. The leading satellite had the z perturbation 
subtracted from the nominal initial conditions, increasing the satellite’s velocity. The 
trailing satellite had the z perturbation added, slowing its velocity. The result was that 
the satellites were drifting apart at about twice the rate they would if only one satellite 
had the perturbation and the other satellite remained unchanged, as was the case for the 
first group (Table 4.4). 

From Table 4.4, it can be seen that for the range of the given perturbations, the 
results appear to be linear. That is, an order of magnitude increase in the perturbations 
results in an order of magnitude change in latitude errors and drift rate values. The drift 
rate between the satellites seems to only depend upon the total distance the satellites are 
displaced from the nominal. If both satellites are in error by the same amount and in the 
same direction, no drift rate will be generated. The error in repeatability can be adjusted 
by correcting the error in position, but the drift rate can be eliminated by either 
correcting the offending satellite, or by causing an equivalent error to occur in the 
repeating satellite. It also seems to be irrelevant whether the perturbations are positive 
or negative as far as the magnitude of the resulting latitude or drift rate errors are 
concerned. A positive change in latitude indicates that the satellite is ahead of its 
desired final position; a negative change, indicates that the satellite is behind its desired 
final position. Though this study was conducted for satellites in a 9 x 9 geopotential 
field, studies presented in Chapter 5 indicate that this table is reliable regardless of the 
size of the static geopotential field being used to determine the motion of the satellites. 
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4.6 Range of Reliability of the Linear Approximation 

Using the nominal conditions for the leading satellite in the 9 x 9 geopotential 
field, an analysis of the linear range of the perturbations is presented. This analysis is 
to ascertain if a large deviation in the radial direction of the initial conditions will alter 
the orbit's ability to remain frozen, and to insure that a perturbation in the initial 
conditions will result in a purely linear change in the drift rate and in the latitude error. 


The nominal initial conditions were given a perturbation in the positive z- 
direction (essentially radially) of a specified amount as indicated in Table 4.5. The 
resulting differences in the final conditions from those of the nominal orbit were 
computed. The differences of each coordinate were squared and summed, and the 
square-root of the result was used as the deviation between the two cases. 


Table 4.5 

Linear Reliability of Changes to the Initial Conditions 


Perturbation in the Initial z-Component 
0.01 meters 
0.1 

1.0 

10.0 

100.0 


Position Deviation After 32 Sidereal Davs 
98.98 meters 
989.42 
9894.31 
98943.1 
987113.5 


1000.0 


8973198. 
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The change in the nominal initial conditions were increased by an order of 
magnitude, starting at one centimeter and finishing at one kilometer. Table 4.5 contains 
the results. As long as the deviations in the final conditions remain proportional to the 
change given to the initial conditions, the error will remain within the linear region. 
With a change in the initial radial distance of 1000 meters, the errors are entering the 
nonlinear region. If the errors in the initial conditions are larger than 1000 meters, then 
the technique derived in Section 4.3 to correct the initial conditions may not be reliable. 

Even with a perturbation as large as one kilometer, the orbit remains frozen, 
since one kilometer is a small percentage change in the semimajor axis. The limitation 
on the maximum perturbation before the orbit becomes non- frozen was not investigated 
since an error as large as one kilometer is unreasonable as far as the mission 
requirements are concerned. The orbits will remain frozen for several orders of 
magnitude beyond the acceptable limits for the error permitted in the satellites' 
positions, so the frozen orbit appears to be a very stable configuration. 

The question of stablity of the frozen orbit due to nonconservative forces was 
not investigated in this study, however, the effects of atmospheric drag were 
investigated in works by Nickerson, et al. [1978] and by McClain [1987]. They 
concluded that atmospheric drag was a devastating problem in maintaining the frozen 
orbits integrity; however, this effect should not concern the GRM satellites since they 
will have drag compensation system on board. 
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4.7 Correction of Velocity 

As described in Section (2.3), the satellites will be initially inserted into a 275 
km altitude orbit [Keating, et al., 1985] and after a series of maneuvers, the satellites 
will descend to the operational altitude of 160 km. Orbit corrections will be required to 
place the spacecraft into the proper mission orbit. Once the satellites are in orbit, an 
instantaneous adjustment in position to eliminate the errors in the initial conditions is 
not be possible. Instead, an adjustment in their velocities that will correct for any orbit 
errors will have to be determined. 

When an orbit is frozen and polar, there are no secular or long period trends in 
any of the orbit elements, as is in the case of GRM. The only gravitational effects that 
appear in the orbital elements for a frozen, polar orbit are short period effects. The 
mean orbital elements are constant, are not influenced by long period effects, and 
excluding the short period trends, the orbit has the behavior of a two-body orbit. Since 
there are no nonconservative forces acting on the proof masses' trajectories, then 
energy is conserved. Energy can be approximated as the two-body energy using mean 
orbital elements: 


E = v 2 /2 - (i/r (4.3) 

where E is energy, v is the orbital velocity and r is the mean orbit radius. Because the 
mean orbital elements for a frozen orbit are constant, then excluding the short period 
effects, energy is constant and the variation in the energy equation is then equal to zero. 
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8E = 0 = v 8v 0 + |i/r 2 8r 0 (4.4) 

Again, since the adjustment will be made near the pole, z can be used instead of 
radius or semimajor axis. With the initial configuration of the satellites being close to 
the pole, and with the satellites traveling mainly in the y-z plane, a change in the z 
position produces a change in the y component of the velocity. Therefore, the change 
in velocity, 8v 0 , will be 8y 0 . From Equation (4.4), a change in the z position, will 
correspond to a change in the velocity as: 

8y 0 = -|V(a 2 v) 8z 0 (4.5) 


Since the orbit is nearly circular, velocity will be approximately equal to the following: 

8y 0 = -M/(a 3 n) 8 zq (4.6) 

The two-body mean motion is equal to (p/a 3 ) 1/2 , which will simplify the change in the 
y component of velocity to: 


8y 0 = -n 8z 0 (4.7) 

This relationship between the change in the initial position to the change in the 
initial velocity was also determined by Colombo [ 1984]. From Equation (4.7), setting 
8z 0 equal to one meter is equivalent to an adjustment in the initial velocity of 
0.001 19636 m/sec in the positive y direction. For verification, a positive change in the 
z direction of one meter was added to the leading satellite. This resulted in a drift 



between the two satellites of about 300 m/day and caused the leading satellite to no 
longer have an exact repeat (0.086° error in latitude). With the adjustment in the initial 
velocity given above, the drift decreased to about 5.6 m/day and the error in latitude for 
the leading satellite was reduced to 0.001°. 

Given that there is a drift between the two satellites and that they do not exactly 
repeat, a method for adjusting the initial position to correct for these discrepancies was 
developed in Section (4.3). Since, during flight, the satellites' positions cannot be 
instantaneously altered, this new procedure will transform the needed positional change 
in the z-coordinate to a correction in the y component of the velocity which, unlike 
positional changes, can be adjusted. 

4.8 An Earlier Determination of Orbital Adjustments 

The technique employed to obtain initial conditions that produce an exact 
groundtrack repeat requires that the mission continue for thirty-two sidereal days. 
Since the initial conditions must repeat after thirty-two days, the desired positions of the 
satellites are known at that future time; the position after 32 days can be compared to the 
initial conditions and the appropriate adjustments made. This procedure does not 
require a priori knowledge of the geopotential, which is convenient since the 
determination and improvement of the geopotential is a primary goal of this mission. 
Waiting for thirty-two days, however, may not be very practical for the actual mission, 
and indeed, once the satellites are operational, only seven days will be available to 
correct any discrepancies in the orbit's iritial conditions [Keating et al., 1986]. 
Consequently, another technique must be used. 
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If a repeating groundtrack could be simulated that would be close to the actual 
groundtrack, then the errors in the satellites' positions could be approximated within the 
first week of the mission instead of waiting for the entire groundtrack repeat cycle. The 
deviations in the actual groundtrack from the designed groundtrack, after only a few 
days into the mission, could be used to determine the proper initial condition 
adjustments that would produce an exact groundtrack repeat for the actual mission. The 
expectation is that, after only a week, the differences between the designed and the 
actual groundtracks will not have grown too large to invalidate this assumption. 

This approach was attempted by comparing a groundtrack generated from a 9 x 
9 geopotential field to a groundtrack from a 36 x 36 geopotential field (both subfields of 
OSU322). Results from this comparison are presented in Figures 4.6a, 4.6b and 4.6c. 
After only seven days, the differences between the two groundtracks, for the leading 
satellite, have grown to 1.6 kilometers (Figure 4.6c). These differences are probably 
generated from the order 16 and order 33 harmonic terms, which are in resonance with 
the GRM satellite altitude, and are not present in the 9 x 9 field. 

The required corrections to the 9 x 9 geopotential’s initial conditions to obtain a 
repeating groundtrack with the 36 x 36 geopotential .are 8 zi = 3.845514 m and 8 z 2 = 
2.611069 m, determined using the full groundtrack cycle for this simulation of 32 
sidereal days. Figures 4.7a and 4.7b show the differences between the actual required 
corrections to the initial conditions, and the corrections determined after the specified 
number of days for each satellite. These two plots, for each satellites' corrections, are 
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nearly identical. These figures have an oscillation about the required solution with a 
bias of approximately a half a meter . 

A mean error offset of 50 cm is too large to completely eliminate the errors in 
the initial conditions, however, these corrections are 85% of the required corrections 
and could be used to reduce some of the errors in longitude and decrease the magnitude 
of the drift rate. The groundtracks of the actual and the modeled trajectories will have 
to be closer to make a proper comparison. To achieve a closer groundtrack, the 
principal resonant terms must be included in the model of the simulated groundtrack. 

A more representative example for a groundtrack comparison is presented using 
the 180 x 180, OSU322 [Rapp, 1981] geopotential field and the 360 x 360, OSU86F 
field [Rapp and Cruz, 1986]. The OSU322 field was used to generate the two 
satellites' reference orbits that repeat after 32 sidereal days to within two kilometers. 
Chapter 5 describes this simulation in more detail. The final conditions from the 
OSU322 simulation indicate a closure for the leading satellite of 1.87 kilometers, and a 
closure of only 111 meters for the trailing satellite (Table 5.2b). The observations, or 
"actual" orbit points, were generated from the OSU86F field which had a groundtrack 
closure of over 154 kilometers each. These two geopotential fields have no common 
harmonic coefficients, and even have different gravity parameters ( IA|i| = 2.x 10 8 
m 3 /sec 2 ), as well as different mean equatorial radii ( lAael = 8 m ). In addition, unlike 
the previous groundtrack comparison example, these two simulations had different 
initial conditions. 


This comparison is more representative of the actual mission since the true 
satellite motion will more closely follow a path modeled by a 180 x 180 field than it 
would the 36 x 36 subfield. In addition, there are no significant resonant terms higher 
than order 180, where the exclusion of resonant terms was a problem with the previous 
comparison. Therefore, the trajectories generated from the entire OSU322 field were 
expected to be comparable to the trajectories generated from the OSU86F field. 

The initial conditions and the results from the simulation generated using the 
OSU322 field are presented in Table 5.2a. The figures of the relative motion for this 
simulation are also presented in Chapter 5 (Figures 5.1a and 5.1b). This simulation 
was used as the design orbit which had a repeating groundtrack. The "actual" trajectory 
was based on the OSU86F field and did not have a repeating groundtrack. Figures 
4.8a and 4.8b illustrate the relative motion for the "actual" trajectory. The initial 
conditions and the results from the OSU86F trajectory are provided in Table 4.6a and 
the initial and final latitudes and longitudes are provided in Table 4.6b. 

From the linear technique described in Section 4.3 calculated using the entire 32 
sidereal days, the actual corrections needed in the initial conditions for the satellites in 
the OSU86F field to cause them to repeat to within the one kilometer were 16.06954 
meters for the leading satellite and 16.62127 meters for the trailing satellite. The 
differences in latitude as a function of time between the repeating, design orbit and the 
observed, actual orbit were computed, from comparisons made every 400 seconds. 
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Table 4.6a 

Initial Conditions for OSU86F Simulation 
Nonrepeating Groundtrack 




Position (mi 


Velocity(m/s) 

Leading Satellite 

r x 

262.16184162 

v x 

-0.048185197 


r y 

-150104.5682242 

v y 

-7816.577574349 


r z 

6515208.810795 

v z 

-179.577052647 

Trailing Satellite 

r x 

262.89992177 

Vx 

-0.0447663774 


r y 

149884.9023112 

v y 

-7816.570937516 


r z 

6515238.3234627 

v z 

179.5005234656 



Final 

Conditions for OSU86F Simulation 
Nonrepeating Groundtrack 

Leading Satellite 

r x 

258.98730736 

v x 0.02556299312 


r y 

-308288.03591398 

v y -7810.9707021185 


r z 

6508774.1727936 

v z -369.3070304526 

Trailing Satellite 

r x 

256.16572671 

Vx -0.0431347331 


r y 

-14138.05881167. 

v y -7819.700419631 


r z 

6516015.8600261 

v z 172.9347354153 
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Table 4.6b 

Latitude and Longitude Values for each Satellite 
from the OSU86F Simulation 


Initial 

Latitude 

Longitude 


Lead Satellite Trailing Satellite 

88.6888° 88.69072° 

169.76012° 349.55956° 


Final after 32 Sidereal Days 

Latitude 87.3059° 89.87647° 

Longitude 169.70819° 170.6999° 


The required changes in the initial conditions in the radial direction are provided in 
Figures 4.9a and 4.9b for each satellite. The true drift rate between the two "actual" 
satellites was 151 m/day, determined from the relative range for the entire 32 days 
(Figure 4.8a). The drift rates used for the daily predictions were computed with data as 
it was accumulated, at the end of each sidereal day the drift rate was updated. These 
drift rates are tabulated in Table 4.7. 
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Table 4.7 

Drift Rate Values for the Given Number of Days 
Davs Drift Rate (m/day) 


1 

2 

3 

4 

5 

6 

7 

8 


296.8 

251.622 

205.94 

180.996 

187.3 

196.29 

192.27 

178.685 


The mean values from Figures 4.9a and 4.9b after eight days were 16.057 
meters for the leading satellite and 16.775 meters for the trailing satellite. These results 
indicate that a repeating groundtrack generated by a 180 x 180 field can sufficiently 
predict the corrections needed for a significantly largo* gravity field in under 8 days. 
By using these mean values from the daily calculations to correct the initial conditions 
of the actual orbits, there will be a resulting groundtrack closure of under one kilometer 
for each satellite. 
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4.9 Summary 

In this chapter, it has been shown that corrections to the satellites' initial 
conditions can be predicted after one complete repeat of the groundtrack. With only 
one week of observations, accurate predictions of the needed corrections can be 
determined if a suitable design orbit is available. During the actual mission, a one week 
set of data can be collected and the calculated corrections to the initial conditions can be 
numerically integrated forward to the current mission time to determine the corrections 
to the satellites' orbits. 

The corrections to the initial conditions are linear if the perturbations are 
generated by a static geopotential field. In Chapter 6, temporal perturbations to the 
satellites' motion are investigated. The errors in the latitude and longitude produced by 
the temporal perturbations may not comply with the linear prediction of latitude and 
longitude errors, or the required radial adjustments that were presented in Table 4.4. 
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Figure 4.1a 
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Figure 4.2b 

Change in y-coordinate for an orbit 
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Figure 4.4a 
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Figure 4.4c 
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Figure 4.5a 

Relative range difference from 300 km . Drift rate was 5.9 m/day 
For 64 sidereal days: Non-frozen orbit 
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Figure 4.5b 
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Figure 4.6c 
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Difference in the actual connection needed and the calculated 
correction versus time for a 36 x 36 geopotential field. 
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Figure 4.9a 

Difference in the actual correction needed and the calculated 
correction versus time for a 360 x 360 geopotential field. 
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Figure 4.9b 

Difference in the actual correction needed and the calculated 
correction versus time for a 360 x 360 geopotential field. 





CHAPTER 5 


ANALYSIS OF GEOPOTENTIAL RESEARCH MISSION SIMULATION 
5.1 Introduction 

The results of a Geopotential Research Mission simulation are presented in this 
chapter. The simulation spans a 32 sidereal day mission lifetime, the time interval 
selected for one complete groundtrack repeat. An analysis of the orbit residual errors 
and a nominal gravity model to reduce the residuals as specified by the mission 
requirements given in Section 2.5 have also been determined. 

Identification of the dominant resonant coefficients contributing to the satellites' 
motion was investigated. Significant effects due to resonance were found to result 
from spherical harmonics of degree and order greater than 36. The effect that the 
resonant terms have on the relative motion and on the repeatability of the groundtracks 
is discussed as well. 

This simulation considered artificial measurement data that can be used to test 
the proposed techniques, discussed earlier in Section 2.3, to recover the Earth's gravity 
field. The measurement data is the relative range-rate between the two satellites. 
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Although only the gravational forces generated by the Earth are included in the 
simulation described in this chapter, the effects from other forces, as well as kinematic 
effects will be discussed in Chapter 6. 

A preliminary simulation of the GRM satellites' motion was performed on the 
CRAY X-MP/48 computer located at Cray Research Incorporated, Mendota Heights, 
Minnesota, in November, 1985. This simulation used Colombo's initial conditions, 
given in Chapter 4. The results from this earlier simulation were described by Schutz, 
et al. [1986]. The closure of the groundtracks in the Mendota Heights simulation were 
2.36 km and 4.59 km for the leading and trailing satellites, respectively, which were 
well within the 10 km closure criteria. It was found that the satellites drifted apart at the 
rate of 93 m/day in that simulation. A new simulation was performed for this study 
which used an improved set of initial conditions with the expectation that an improved 
groundtrack closure and a reduction in the drift rate between the two satellites would 
occur. This new simulation is described in the following sections. 

5.2 Description of Simulation 

The numerical computations for the new simulation were performed on the 
CRAY X-MP/24, located at the University of Texas System Center for High 
Performance Computing. The amount of computer time required for the simulation 
was approximately 5 hours and 40 minutes for the 32 sidereal day simulated mission. 

The numerical technique for integrating the satellites equations of motion was 
the Encke method, described by Roy [1978] and Lundberg [1985], which used a 
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reference orbit based on a secularly precessing ellipse with an analytical representation. 
The difference between the analytical reference orbit and the true orbit is referred to as 
the "Encke vector", which was integrated in place of the actual satellite state. A primary 
advantage of this method is that it reduces round-off errors associated with numerical 
integration [Lundberg, 1985]. The reference orbit characteristics are provided in Table 
5.1, and were selected to produce a small magnitude Encke vector without secular 
trends. 


Table 5.1 

Secularly Precessing Reference Orbit for Encke Method 


Lead Satellite 
a = 6523600.811305 m 
e = 0. 
i = 90° 

Cl = 90° 
co = 0° 

M= 1.59331 106462 rad 
Cl = 0 rad/day 
(0 = 0 rad/day 

M = 0.001 1963632130 rad/day 


Trailing Satellite 
a = 6523599.627289 m 
e = 0. 
i = 90° 

Cl = 90° 
co = 0° 

M= 1.547312542033 rad 
Cl = 0 rad/day 
<o = 0 rad/day 

= 0.001 19636336526 rad/day 
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The integrator is a class 2, for second order ordinary differential equations, 
fixed-mesh, multistep algorithm, of order 10 as described by Lundberg [1985]. The 
integration step size was five seconds, chosen to be small enough to detect the highest 
degree harmonics used in the simulation [ Schutz , et al., 1986]. 

The 180 x 180 OSU322 gravity field described in Section 1.2 was modified by 
including terms out to degree 300 and lower order harmonics to 10. The force model 
used the gravity parameter of 3.9860064 x 10 km /sec and the mean equatorial radius 
of 6378.145 km. The epoch time was chosen to be 2445700.5; midnight, January 1, 
1984, to be consistent with the previous simulations. The sidereal hour angle was 
100.1 135613° and the Earth's rotation rate was held constant at 7.2921 1585531 x 10" 5 
rad/sec. 

For consistency with the previous work in this report and the adopted models, a 
new set of initial conditions was calculated for the simulation. These new initial 
conditions were based on the updated values of the gravity parameter and the mean 
equatorial radius, as well as the 36 x 36 OSU322 subfield. These initial conditions, 
calculated by the method described in Section 4.3, were expected to lead to a smaller 
groundtrack closure than resulted from the initial conditions derived by Colombo. For 
the groundtrack to repeat to within 10 km, the error in the geodetic latitude must be less 
than 0.1°. 

The instantaneous relative range was calculated by subtracting the two satellites' 
instantaneous states at each time point. The relative range vector is: 



101 


P=r 1 -r 2 (5.1) 

where r is the position vector of the leading satellite and i* 2 is the position vector of the 
trailing satellite. The relative range-rate, p, is calculated by: 

p = p -p (5.2) 

P 

Figures 5.1a and 5.1b show the relative motion between the two satellites from 
the simulation. Figure 5.1a illustrates the relative range difference which was 
determined by subtracting the actual distance between the two satellites from their 
desired average separation distance of 300 km. A nonlinear trend was present, which 
was due mainly to resonant terms of the order 82 terms, though other resonant terms 
have a contribution. Because the initial conditions were created to be an exact repeat 
with no drift for a 36 x 36 gravity field, a secular trend exists in the relative range 
measurements when the higher harmonics were included. For this same reason, the 
satellites did not repeat exactly. The initial and final conditions for both satellites are 
provided in Table 5.2a. The satellites had a drift rate of approximately 41 m/day. The 
error in the groundtrack repeat was less than two kilometers for the leading satellite and 
approximately 100 meters for the trailing satellite; well within the closure criteria 
specified earlier. The initial and final geodetic latitude and longitude after 32 sidereal 
days are presented in Table 5.2b. The drift rate between the satellites and the latitude 
errors can be eliminated by an adjustment in the initial conditions (as described in 
Chapter 4), but the periodic trend due to the resonant terms cannot be removed. 
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Table 5.2a 

Initial Conditions for OSU322 Simulation 
Repeating Groundtracks 




Position (m) 


Velocitvfm/sl 

Leading Satellite 

r x 

253.7524 

v x 

0.000000 


r y 

-150000.0000 

v y 

-7816.570937552 


r z 

6515237.772962 

v z 

-179.497721338 

Trailing Satellite 

r x 

253.7524 

v x 

0.000000 


r y 

150000.0000 

v y 

-7816.570937516 


r z 

6515240.934532 

v z 

179.5005234656 


Final Conditions for OSU322 Simulation 
Repeating Groundtracks 


Leading Satellite 

r x 

255.2931288 

v x 

0.0799469041 


r y 

-147951.4011902 

v y 

-7818.18594707 


r z 

6513972.4970525 

v z 

-177.5348494778 

Trailing Satellite 

r x 

253.5177113 

V x 

0.0214660832 


r y 

149988.6625145 

v y 

-7818.19069170 


r z 

6513896.6701621 

V Z 

179.0361331515 
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Table 5.2b 

Latitude and Longitude Values for Each Satellite 


Initial 

Lead Satellite 

Trailing Satellite 

Latitude 

88.6888° 

88.69072° 

Longitude 

169.76012° 

349.55956° 

Final after 32 Sidereal Davs 


Latitude 

87.707° 

88.689° 

Longitude 

169.759° 

349.563° 


5.3 Investigation of the General Behavior 

The simulation specification that the satellites' orbit complete 525 revolutions in 
32 sidereal days yields a period for the satellites of approximately 88 minutes. In one 
sidereal day, the satellites will have completed about 16.40624 revolutions and the 
orbits will be in resonance with any harmonic terms whose periods are commensurate 
with integer multiples of the satellites’ daily revolutions [Kaula, 1966]. 

From Equations (3.2), the periods of the resonant terms can be determined. 
The denominator of these equations is a function of the frequency associated with a 
particular harmonic coefficient. The frequency equation for a particular perturbation is: 
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'F = (/-2p)d> + (l-2p+q)M + m(Q-€>) 


(5.3) 


where the parameters are defined earlier in Chapter 3 [Kaula, 1966]. If this 
denominator term approaches zero, the orbit is in resonance with the harmonic term of 
degree /, and order m. For this particular mission, the time rate of change for the 
argument of perigee, d), is equal to zero because the orbit is frozen; and the nodal rate, 
Q, is zero because the orbit is polar. The only remaining terms are the ones associated 
with time rate of change of mean anomaly and the rotation rate of the Earth. The closer 
the integer multiples of M are to the values for the integer multiples of sidereal days, the 
deeper the resonance associated with the particular harmonic term becomes, and the 
more significant that term will be in the satellite's motion. 

Table 5.3 illustrates the order of the harmonic terms whose periods are in 
resonance with the GRM satellites' periods. The terms that generate the deepest 
resonance are under order 180. The resonant terms, in order of the largest individual 
effect, are orders 82, 33, 49, 164, 115, and 16. The order 82 terms were dominant, 
with an amplitude of over 800 meters and have a period of approximately 32 days. The 
effects generated by each of the dominant resonant terms are investigated separately. 
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Table 5.3 

Resonant Terms Plus Side Bands 
4* = (/-2p)d) + (l-2p+q)M + m(Q-©) 
n (16.4062414) - m = 1/D 
n = l-2p+q 


a. 

tn 

m 

D CDavs) 

1 

16 

0.4062414 

2.46 

1 

17 

0.5937586 

1.684 

2 

32 

0.8124828 

1.230795 

2 

33 

0.1875172 

5.33284 

3 

49 

0.2187242 

4.57196 

3 

50 

0.7812758 

1.279957 

4 

65 

0.6249656 

1.6 

4 

66 

0.375035 

2.666666 

5 

82 

0.031207 

32.04409 

5 

83 

0.968793 

1.032212 

6 

98 

0.43744 

2.286027 

6 

99 

0.5625516 

1.777614 

7 

114 

0.8436898 

1.185269 

7 

115 

0.1563102 

6.397535 

8 

131 

0.2499 

4.00 

8 

132 

0.7500688 

1.333211 

9 

147 

0.6561726 

1.523989 

10 

164 

0.062414 

16.022 

10 

165 

0.937586 

1.066568 
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With the same initial conditions used for the simulation, the 36 x 36 OSU322 
subfield was the baseline geopotential field and individual order resonant terms were 
added to observe their effects. Only the five dominant terms associated with harmonics 
of resonant order were investigated. The relative motions plots of the results are 
included, and Table 5.4 contains the results of the effects the resonant terms had on the 
final longitude and latitude values. The zonal harmonics to 300 were also included in 
this investigation to insure that no long periodic effects due to odd zonal harmonics 
remain. 


Table 5.4 

Effects on Final Latitude and Longitude due to Resonance 


Order 

82 

33 

49 

115 

16 

164 

zonals to 300 


Satellite 1 
8(lat°/long°) 
0.172/0.01 
-0.123 /-0.009 
0.057 / 0.003 
0.001 / 0.0 
0.007 / 0.0 
0.003 / 0.0 
-0.029 /-.003 


Satellite 2 
8(lat°/long°) 
-0.146/0.01 
0.134 /-0.008 
-0.052 /-0.003 
-0.004 / 0.0 
-0.005 / 0.0 
-0.006 / 0.0 
0.031 /-0.002 


Drift rate 

(m/day) 

-116.44135 

-36.07876 

-14.85071 

15.09051 

-1.828319 

14.789918 

-0.0243366 


Having a period of 32.044 days and an amplitude of 806 meters, the harmonic 
terms of order 82 caused the deepest resonant effect. With a 36 x 36 plus order 82 
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geopotential field, the satellites had latitude errors of 0.172° for the leading satellite and 
0.146° for the trailing satellite. A drift between the satellites of -116.44 m/day was 
produced. The relative motion is illustrated in Figures 5.2a and 5.2b. The 32 day 
period cannot be seen clearly in the relative range because of the large secular trend. 
With a slight adjustment in the initial conditions, however, the relative range history can 
be adjusted to have the characteristics similar to the full field's relative range (Figure 
5.2c). This similarity demonstrates that the order 82 resonance was the major 
contributor to the nonlinear trend in the relative range shown in Figure 5.1a. 

The second most dominant resonant harmonics were the order 33 terms 
(Figures 5.3a and 5.3b). The period associated with this order was only 5.33 days 
but, they produced an amplitude of 83 meters in the relative range. These harmonics 
caused the satellites to lag behind a repeat groundtrack by 0.123° in the leading satellite 
and 0. 134° in the trailing satellite. The satellites drift apart at the rate of 36.08 m/day. 

The order 49 terms, with a period of 4.57 days, were considered the next most 
dominant (Figure 5.4a and 5.4b). Even though they had an amplitude of only 19 
meters, they caused a deviation in the groundtracks of 0.057° and 0.052°; which was 
the third largest change in latitude. The drift rate caused by adding this order to the 36 
x 36 field was -14.85 m/day. 

The third highest amplitude of 40.65 meters was due to the resonant terms of 
order 164 terms (Figure 5.5a and 5.5b). These terms have a period of 16.022 days, 
half the period of the order 82 terms. However, the effect on the the groundtrack was 
only 0.003° in latitude for the leading satellite and 0.006° for the trailing satellite. The 
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drift rate from these terms cause the satellites to approach each other at the rate of 14.85 
m/day. Although these terms had little effect on the orbits repeatability, the magnitude 
of the amplitude of their oscillation was relatively large, thereby causing these terms to 
be significant, even though they are of high degree. 

The order 115 terms had a larger effect on the satellites motion than the order 16 
terms, but both terms contributed very little compared to the effects of orders 82, 33, 
49, and 164. The order 16 terms were included because they are within the nominal 36 
x 36 field and do not add to the overall size of the geopotential fields when included, 
and the order 115 terms were included because they had a larger effect than the order 
16. The order 16 terms have a period of 2.46 days, with an amplitude of 10.12 meters 
(Figure 5.6a and 5.6b). The order 1 15 terms have a 6.4 day period and an amplitude 
of 17.64 meters (Figure 5.7a and 5.7b). 

Finally, the effect of the zonal harmonics up to degree 300 were investigated to 
insure that no secular trend was generated by the odd zonal harmonics. The satellites 
geodetic latitude disagrees with the repeating path by -0.029° and -0.031° for each 
satellite. This equates to an error of about three kilometers on the Earth's surface, and a 
drift rate of less than -0.024 m/day was incurred between the satellites. The odd zonal 
harmonics produced a long period of 79 days, thus an effect with this period might 
appear to be secular in the in a 32 day time interval. However, with a drift rate of only 
0.024 m/day, the amplitude of this long period oscillation will be quite s mall , yielding a 
maximum possible amplitude of only two meters. The conclusion can be made that the 
long periodic effect due to the odd zonal harmonics have been eliminated by selecting a 
frozen orbit that was based on zonal harmonics to J9 only. 
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5.4 Reduction of Residuals 

The simulated ephemerides for one of satellite described in Section 5.2 was 
taken as a set of observations for the program UTOPIA, which performed a least 
squares fit of the observation data. The observations were the set of inertial (J2000), 
geocentric values of the position vector of the leading satellite, Y,-, where i is the time of 
the observation. The calculated set of observations are G( X, ,t, ), where X,*, is the 
nominal state of the satellite as it travels along the path of the orbit determined by the 
specified geopotential model [Tap ley, 1972]. The difference between the observations 
for the leading satellite, provided by the simulation, and the calculated set of 
observations of the leading satellite are the observation residuals, r, : 

r,- = Y,-G(X i *,t < - ) (5.4) 

To reduce the magnitude of the residuals, either the model can be improved by 
altering the values assigned to the gravity coefficients, expanding the model by 
including more terms that contribute to the satellite's motion, and/or by changing the 
initial conditions. Since the simulation that generated the set of observations used only 
the forces generated by the Earth's geopotential field, these will be the only forces 
included in the nominal model. 

The classical orbit determination technique of least squares was used to 
determine a nominal geopotential field that reduced the residuals to within the mission 
specifications. The specifications, given in Section 3.4, are 100 meters in the radial 
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direction, 300 meters in the transverse direction, and 300 meters in the normal direction 
[Keating, et al., 1986]. These are the largest permissible errors in the satellites' 
positions that will allow the recovery of the Earth's geopotential field to the desired 
resolution. In addition, 99.98% (or 3a) of the residuals must be within these 
specifications. It is desired to use the smallest geopotential field for the nominal model 
that will cause the residuals to be under the mission specifications in order to minimize 
the computational effort. The nominal geopotential model was designed to reduce the 
residuals to a level that a more efficient recovery technique could be used to further 
refine the geopotential model. 

To reduce the residuals, the initial conditions can be estimated by performing a 
least squares fit with the observation data None of the harmonic coefficients were 
estimated, and no resonant terms were included in the baseline model. The radial 
component of the residuals varied between ±200 meters and the normal component 
varied between ±100 meters. The transverse component of the residuals varied between 
±4500 meters, indicating that the estimation of the initial conditions alone will be 
inadequate for a 36 x 36 geopotential field (Figure 5.8a, 5.8b and 5.8c). From the 
Figure 5.8b, a 32 day period can clearly be seen, along with an approximately 5.5 day 
period superimposed. These two periods were caused by resonant coefficients of order 
82 and order 33. Since the residuals were not within the mission specifications, the 
coefficients in the gravity model must also be estimated. 

It is desired to use the minimum number of coefficients above the 36 x 36 field 
as possible. The resonant coefficients, in order of their greatest effect that were 



candidates to be estimated were 82, 33, 49, 164, and 16. Only the first two pairs of the 
harmonic coefficients were used, which when "tuned" along with estimating the 
position and velocity, produce the best fit to the observation file generated by the full 
OSU322 gravity field simulation. With the simulation and the estimated model having 
the same OSU322 36 x 36 subfield, the importance of the resonant terms can be 
illustrated. The results from estimating the resonant terms are presented in Table 5.5. 
This table includes the root mean squared values (RMS), as well as the maximum 
magnitudes the residuals obtained. 

The 36 x 36 nominal model was expanded to include the estimation of the first 
two pairs of order 82: Cs2,82 and S82.82, and Cs3,82 and Ss3,82 (Figure 5.9a, 5.9b 
and 5.9c). In the estimation of these four coefficients, not only will they need to 
account for the differences in the two models used, but they will also need to absorb the 
effects of the omitted order 82 coefficients. An a priori covariance was assigned to the 
estimated harmonics using Kaula's rule: 


c= ±\0' 5 

l 

where / is the degree of the harmonic and a is the standard deviation [Kaula, 1966], 
The residuals have been reduced to be within 120, 940, and 70 meters in the radial, 
transverse, and normal directions. Thus, the addition of the order 82 terms in the 
estimation has improved the residuals by a factor of four in the transverse direction, but 
the residuals were still too large, consequently, the order 33 terms were included. 
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Table 5.5 

Reduction of residuals 
OSU322 36 x 36 Geopotential field 


Order 


RMS (ml 

Lareest ( m ) 


R 

32.0 

117.137 

82 

T 

438.12 

940.62 


N 

25.34 

70.73 


R 

28.68 

106.93 

33 

T 

155.31 

444.86 


N 

17.408 

52.79 


R 

27.82 

103.98 

49 

T 

101.42 

314.79 


N 

15.76 

50.755 


R 

27.822 

104.07 

16 

T 

101.25 

309.31 


N 

15.84 

51.97 


R 

27.19 

107.18 

164 

T 

75.63 

312.60 


N 

15.82 

52.36 


The first two pairs of geopotential coefficients from order 33 were estimated 
along with the first two pairs of the order 82 terms. The residuals were reduced by a 
factor of two from estimating order 82 alone (Table 5.5), and for this model, were 
within 107 meters in radial, 450 meters in transverse, and 53 meters in the normal 
directions. Figures 5.10a, 5.10b and 5.10c are the radial, transverse and normal 
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residuals versus time. From Figures 5.10b, that illustrate the transverse residuals, a 16 
day period can be seen, caused by the order 164 terms. Superimposed on the larger 
amplitude 16 day period is a smaller amplitude oscillation from the order 49 terms. 
Since the order 49 of lower order than 164, then to save computation time and to 
further reduce the residuals, the order 49 coefficient pairs were estimated in preference 
to the order 164 terms. 

The estimation of the first two coefficient pairs of order 49, as well as the 
previous harmonic terms, reduced the residuals to within 104 m, 315 m, and 51 m in 
the radial, transverse, and normal directions (Figure 5.1 la, 5.1 lb and 5.1 lc). The 16 
day period is clearly seen in the transverse residuals (Figure 5.11b), but since the 
residuals were close to the mission specifications, estimation of the terms as high as the 
order 164 terms in the model may be unnecessary. Instead, the first two pairs of order 
16 were estimated, because the order 16 terms are within the nominal 36 x 36. 

Estimating the first two pair of these low order resonant coefficients decreases 
the residuals closer to the 3a band required for the geopotential recovery (Figure 5.12 

a, 5.12b and 5.12c). The maximum amplitudes of these residuals were 104 m, 310 m 
and 52 m. The mean values from the residuals are provided in Table 5.5 which indicate 
that on the average, the inclusion of the order 16 terms in the model helps, but to a very 
small degree. 

The order 164 terms, though a deeper resonant order, were the last coefficients 
to be included in the nominal gravity model. Though the maximum amplitudes 
increased somewhat for the residuals, the mean values decreased, and for the transverse 
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residuals, they were decreased by a fourth (Figure 5.13a, 5.13b and 5.13c). This 
indicates that the order 164 terms may need to be included when the baseline 36 x 36 
geopotential field is other than the OSU322 gravity field model. 

Since, in the actual mission, the true geopotential field will be unknown, it 
could be unrealistic to employ the same subfield that generated the simulation to reduce 
the residuals. Therefore, the OSU322's 36 x 36 gravity field was replaced with the 
Goddard Earth Model GEM 1 OB in order to have a perturbation in the baseline gravity 
field. The same gravity parameter of 3.9860064 x 10 5 km 3 /sec 2 and the mean 
equatorial radius of 6378.154 km were used for this analysis. By estimating the same 
gravity harmonics that were estimated using the OSU322 field, the residuals could only 
be reduced to within 500 meters in the radial direction, 1100 meters in the transverse 
direction, and 64 meters in the normal direction (Figure 5.14a, 5.14b and 5.14c). 

It appears that there is too much disparity between the two geopotential fields to 
only estimate the same five pairs of coefficients. To adequately reduce the residuals, 
additional harmonic terms within the 36 x 36 field must be included. The even zonal, 
J2, was then estimated which caused a dramatic decrease in the residuals, but the radial 
term had a rather large bias in the negative direction (Figure 5.15a, 5.15b and 5.15c). 
With the inclusion of the odd zonal, J3, the bias was eliminated, and it further reduced 
the residuals to within the mission specifications (Figure 5.16a, 5.16b and 5.16c). The 
complete set of coefficients that were estimated are: (82,82), (83,82), (49,49), 
(50,49), (33,33), (34,33), (16,16), (17,16), (164,164), (165,164), J 2 , and J 3 . It is 
important to note that the estimated coefficients are regarded as parameters to adjust in 
order to produce a nominal orbit with residual errors less than the mission 
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specifications. These parameters do not represent an improved or even a geophysically 
meaningful set of coefficients. 

From the residuals in the transverse direction, a period of approximately 1.68 
days can be detected (Figure 5.16b). This period was generated by the order 17 terms, 
a side band of the deeper resonant terms, of order 16. Each of the deeper resonant 
coefficients have side bands associated with them that are in weak resonance with the 
satellite's orbital period. For a different simulation, these weaker resonant terms may 
need to be included in the model in order to properly reduce the residuals to the level of 
the mission specifications. 

5.5 Summary 

The conclusions drawn from the experiments described in this chapter indicate 
that there are significant gravity effects due to terms of degree and order greater than 
36. To meet the mission specifications for a nominal geopotential model, the entire 36 
x 36 field does not need to be estimated, but selected higher order terms will require 
adjustments. If the mission specifications become more constrained, then more terms 
may need to be estimated in order to create an adequate nominal gravity model. For the 
error model considered for this study, the dominant resonant terms and the first two 
zonal harmonics are the minimum number of gravity terms needed to establish a 
nominal geopotential model. 

Only the first two pairs of each of the coefficients were estimated. Pairs of 
coefficients were estimated because the odd degree terms have a different frequency 
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than the even degree terms for a particular order. If an attempt is made to estimate more 
pairs of coefficients of the same order then the individual harmonic coefficients may not 
be separable. Instead of estimating more terms of the same order to reduce the 
residuals, another order of coefficients should be included in the model or the a priori 
value for the coefficients could be increased. 

The harmonic coefficients were estimated for the two satellites separately, both 
coming well within the limits given the residuals. Using the leading satellite's 
geopotential field estimated specifically for that satellite's simulated trajectory, only the 
position and velocity of trailing satellite were estimated. The residuals were within 3 o 
specifications of 100 m, 300 m, and 300 m for the radial, transverse, and normal 
components of the residuals. This result indicates that only one gravity field needs to 
be determined and that this field will be suitable for both satellites. 
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Figure 5.1a 

Relative range difference for simulation 
OSU322 gravity field. 32 sidereal days. 
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Figure 5.1b 

Relative range-rate for simulation 
OSU322 gravity field. 32 sidereal days. 
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GEOPOTENTIAL FIELD 0SU322 36 X 36 
RESONANCE TERMS OF ORDER 82 INCLUDED 
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Figure 5.2b 

Relative range-rate 
Gravity field: 36 x 36 plus order 82 
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Figure 5.3b 

Relative range-rate 
Gravity field: 36 x 36 plus order 33 
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GEOPOTENTIAL FIELD OSU322 36 X 36 
RESONANCE TERMS OF ORDER 49 INCLUDED 
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Figure 5.4a 

Relative range difference 
Gravity field: 36 x 36 plus order 49 
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Figure 5.4b 

Relative range-rate 
Gravity field: 36 x 36 plus order 49 
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Figure 5.5a 

Relative range difference 
Gravity field: 36 x 36 plus order 164 
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Figure 5.5b 
Relative range-rate 

Gravity field: 36 x 36 plus order 164 
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GEOPOTENTIAL FIELD OSU322 36 X 36 
RESONANCE TERMS OF ORDER 16 INCLUDED 
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Figure 5.6a 

Relative range difference 
Gravity field: 36 x 36 plus order 16 
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Figure 5.6b 

Relative range-rate 
Gravity field: 36 x 36 plus order 16 
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Figure 5.7a 

Relative range difference 
Gravity field: 36 x 36 plus order 1 15 
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Figure 5.7b 
Relative range-rate 

Gravity field: 36 x 36 plus order 1 15 
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Figure 5.8b 
Transverse difference 

Residuals between simulation and the 36 x 36 gravity field 
from OSU322. No resonant terms were estimated. 
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RESIDUALS BETWEEN FULL FIELD AND 0SU322 36X36 
RESONANCE ORDER 82 ESTIMATED 
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Figure 5.9b 
Transverse difference 

Residuals between simulation and the 36 x 36 gravity field 
from OSU322. Resonant order 82 was estimated. 
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RESIDUALS BETWEEN FULL FIELD AND OSU322 36X36 
RESONANCE ORDERS 82 AND 33 ESTIMATED 
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Figure 5.10b 
Transverse difference 

Residuals between simulation and the 36 x 36 gravity field 
from OSU322. Resonant orders 82 and 33 were estimated. 
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RESIDUALS BETWEEN FULL FIELD AND OSU322 36X36 
RESONANCE ORDERS 82 AND 33 ESTIMATED 
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Figure 5.10c 
Normal difference 

Residuals between simulation and the 36 x 36 gravity field 
from OSU322. Resonant orders 82 and 33 were estimated. 
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Figure 5.11c 
Normal difference 

Residuals between simulation and the 36 x 36 OSU322 gravity field 
Resonant orders 82, 33 and 49 were estimated. 
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Figure 5.13b 
Transverse difference 

Residuals between simulation and the 36 x 36 OSU322 gravity field 
Resonant orders 82, 33, 49, 16 and 164 were estimated. 
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RESIDUALS BET\AEEN FULL FIELD AND OSU322 36X36 
RESONANCE ORDERS 82 33 49 16 AND 164 ESTIMATED 



Figure 5.13c 
Normal difference 

Residuals between simulation and the 36 x 36 OSU322 gravity field 
Resonant orders 82, 33, 49, 16 and 164 were estimated. 
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Figure 5.14a 
Radial difference 

Residuals between simulation and the 36 x 36 GEM10B gravity field 
Resonant orders 82, 33, 49, 16 and 164 were estimated. 
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RESIDUALS BET\^EN FULL FIELD AND GM10B 36X36 
NO ZONALS ESTIMATED 
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Figure 5.14b 
Transverse difference 

Residuals between simulation and the 36 x 36 GEM10B gravity field 
Resonant orders 82, 33, 49, 16 and 164 were estimated. 
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RESIDUALS BET\^EN FULL FIELD AND GM10B 36X36 
J2 AND RESONANCE TERMS ESTIMATED 
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Figure 5.15a 
Radial difference 

Residuals between simulation and the 36 x 36 GEM10B gravity field 
All resonant terms were estimated plus J 2 . 
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CHAPTER 6 


THE EFFECTS OF TEMPORAL PERTURBATIONS 
6.1 Introduction 

The previous chapters described the effects of only the Earth's static gravity 
field on the motion of the GRM satellites, with no other perturbations considered. The 
perturbation study presented in this chapter is only concerned with changes in the orbits 
of the proof masses, therefore, only conservative forces are examined. The effects of 
atmospheric drag, Earth albedo, and solar radiation pressure are excluded by the 
assumption in this investigation that the disturbance compensation mechanism removes 
all nonconservative force effects. 

A comparison was made between the nominal case and the nominal case plus 
the effect due to a specified perturbation. The nominal case was defined in Chapter 4 as 
a groundtrack that repeated exactly, generated with a 9 x 9 OSU322 geopotential field; 
furthermore, there was no drift between the satellites, and the orbits were frozen. The 
difference between the two cases' relative range-rate was used to determine the 
magnitude of the particular perturbation in terms of its affect on the measurement 
signal, and the resultant errors incurred if these effects are omitted from the model. The 
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perturbations in the satellites' motion that were examined are: precession, nutation, 
polar motion, solid Earth tides, ocean tides, lunar, solar, planetary, relativity and the 
Moon's effect on the Earth's oblateness. 

The results from each of these effects will vary depending on the initial 
geometry, or on the epoch time used. Under certain conditions, the perturbations could 
contribute less than the ±1 p.m/sec precision level of the relative range-rate 
measurements (Section 2.4). If a secular or periodic change in the range-rate 
differences exists that results in the magnitude becoming greater than ±1 (im/sec in six 
months time, then that effect must be accounted for in the modeling. 

The software used to make the comparison of the individual effects was 
UTOPIA. The various models used in UTOPIA to generate the effects are defined in 
this chapter. 

6.2 The Effects of Precession, Nutation and Polar Motion 

Precession is caused by the gravitational attraction of the Sun and the Moon 
and to a lesser extent, the planets on the oblate Earth. This gravitational attraction 
causes the Earth's pole to have a westward precession with a period of about 26,000 
years [/toy, 1978]. Precession is the motion of the Earth's mean pole around the 
inertial Z axis, or ecliptic pole [Cappellari, Velez, and Fuchs, 1976]. 

Nutation is caused by the inclination of Moon's orbit to the ecliptic plane and 
the lunar gravitational interaction with the Earth's oblateness. Nutation is periodic 
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rather than secular, and oscillates about the precessional path with a dominant period of 
about 18.6 years and an amplitude of about nine arcseconds (Figure 6.1) [Cappellari, et 
al., 1976]. 

The standard epoch of the fundamental astronomical coordinate system used in 
UTOPIA was Julian Date 241545.0 (January 1, 2000, 12 h ), referred to as Epoch 
J2000.0. The rectangular, inertial, geocentric coordinate system defined by the mean 
equator and equinox of J2000.0 has an X and Y axes located in the Earth's equatorial 
plane. The X-axis is along the vernal equinox of J2000.0 and Z-axis is perpendicular 
to the equatorial plane. This coordinate system is fixed in space, is not influenced by 
precession, nutation or polar motion, and will be referred to as the Mean-2000 
coordinate system. 

The mean-of-date coordinate system is defined by the addition of precession to 
the inertial coordinate system [Cappellari, et al., 1976]. The initial location of this 
coordinate system in space will depend on the epoch date used; for the study reported in 
this chapter, the epoch date was January 1, 1984. The X-axis is directed towards the 
vernal equinox of the epoch date. Because of precession, the equatorial plane will be 
slightly different than the Mean-2000 equatorial plane, therefore, the Z-axis is different 
as well. 

With the inclusion of nutation, the true-of-date coordinate system is defined 
[Cappellari, etal., 1976]. The slight changes in the mean-of-date and the true-of-date 
from the inertial system will effectively shift the orientation of the geopotential field in 
space, and therefore, change the forces acting on the satellites. As with the mean-of- 
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date system, the equatorial plane has been displaced due to precession and nutation. 
When precession and nutation are not considered, then the Mean-2000, mean-of-date, 
and true-of-date systems are all equivalent coordinate systems. 

To investigate the kinematic effects on the satellites' relative motions, the initial 
conditions were formulated in the mean-of-date, or true-of-date (when nutation was 
included), coordinate systems. This choice keeps the satellites in the same relative 
positions to the north pole on the selected epoch date for which they were originally 
derived, thereby producing a better groundtrack repeat than would occur if the initial 
conditions were placed in the Mean-2000 coordinate system. 

The processional effects on the satellites' relative range-rate are shown in Figure 
6.2a. Initially, the difference in the relative range-rate is zero. As time progresses, the 
mean-of-date system moves with respect to the Mean-2000 system, and the changes in 
the force field due to the spatial change in the geopotential orientation begin to 
accumulate. The effects due to the inclusion of precession in the kinematic model 
increase to a magnitude of 35 pm/sec in the first 32 sidereal days. Since this is larger 
than the 1 pm/sec requirement, precession will clearly need to be included in the model. 

Adding nutation to the kinematic model requires the initial conditions to be 
placed in the true-of-date coordinates. Nutation increases the difference in relative 
range-rate to a magnitude of 600 pm/sec, but with an apparent decrease in amplitude 
over the 32 day period (Figure 6.2b). Because nutation oscillates about the 
precessional path, then only when the oscillation intersects the precessional path will 
the true-of-date and mean-of-date systems be coincident. Since the epoch time was not 
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chosen for this to be the case, there is a sudden change in magnitude for the difference 
in relative range-rate when nutation is included in the kinematic model. 

Polar motion, which is the motion of angular velocity vector of the Earth 
relative to the body fixed z-axis, was also investigated and was included along with the 
effects of precession and nutation. There was no discernible effect due to the inclusion 
of polar motion to the model. 

6.3 The Effects of Solid Earth Tides 

The Earth is not rigid, and will deform because of the gravitational attractions of 
the Sun and the Moon. The amplitude of the solid surface deformation can be as high 
as 1/3 meters [ Siry , 1973]. The Earth's geopotential will be altered due to the 
deformation and can be expressed as: 


AU(r) = E k n (ae/r) 2n+1 V n (r) 
n=2 

where AU(r) is the change in the geopotential field at position r, k n is the Love numbers 
of degree n, a e is the mean equatorial radius, and V n is the potential due to the solid 
Earth tides [Shum, 1982]. The Love parameters are an indication of the Earth's 
deformation properties. The solid Earth tide model employed by UTOPIA uses the 
equations for the changes in the geopotential coefficients due to the tidal effects 
provided by the MERIT Standards [1983]. The MERIT Standards uses Wahr's theory 
to model the Earth tides, which uses the 1066a Earth model of Gilbert and Dziewonski 
[1975]. 
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The errors in the relative range-rate due to the solid Earth tides grow to an 
amplitude of 240 pm/s in 32 sidereal days (Figure 6.3). Consequently, the magnitude 
of the differences in the relative range-rate indicate that the solid Earth tides will have a 
significant effect on the satellites' motion and must be modeled. 

6.4 Effects of the Ocean Tides 

The models for the solid Earth and the ocean tides are treated separately. The 
lunar and solar gravitational attraction on the oceans considers the Sun and the Moon to 
be point masses [Torge, 1980]. The ocean tide model used in UTOPIA is based on the 
Schwiderski tide model [1980], which contains a 1° x 1° grid of the amplitudes and 
phase angles for nine of the ocean tide constituents: M 2 , S 2 , N 2 , K 2 , Ki, Oi, Pi, Q, 
and Mf [Eanes, Schutz and Tapley, 1983]. The ocean tide effects are computed by 
summing over the constituents listed above. Values for the amplitudes and phases for 
each constituent are tabulated in the MERIT Standards [1983]. The Schwiderski tide 
model along with modifications to account for the effects of the atmosphere, and the 
expressions for the variations in the geopotential coefficients used by UTOPIA are 
provided by Eanes, et al. [1983]. 

The potential due to the ocean tides contains a 14 day period due to the Moon 
and a 180 day period due to the Sun. The 14 day period can be seen in the relative 
range-rate plot of the ocean tidal effects (Figure 6.4). The tidal potential also contains 
short periods of diurnal and semidiurnal lengths [Torge, 1980]. The change in the 
ocean’s mass distribution generates an effect on the two satellites' relative range-rate. 
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the amplitude of which grows to ±60 fi.m/s in 32 sidereal days. 

A study of the ocean surface variability effects, such as eddies, on the satellites' 
relative range-rate was made by McNamee [1986]. This study concluded that the ocean 
currents will affect the satellites' motion to a maximum value of 20 |im/sec which is 
above the ±1 (im/sec level, and therefore, will also need to be considered. 

6.5 Planetary Effects 

The gravitational effect of the all the planets (Figure 6.5a) on the GRM satellite 
motion was considered. Of all the planets, Jupiter caused the largest change in the 
relative range-rate (Figure 6.5b). The planets were assumed to be point masses in the 
force model. The perturbation to the two-body force for N-bodies is given as: 

F = - £ GMi ( A i/Ai 3 - r-Jx? ) (6.1) 

i=l 

where F is the force on the satellite due to the body Mj, M is the mass of the Earth, A i 
is the vector from the perturbing body to the satellite, and rj is the vector from the 
perturbing body to the Earth. The UTOPIA software uses a planetary ephermerides 
that provides the values for A \ [ Shum , 1982]. The ephermerides for the planets, as 
well as the Sun and the Moon, are from the Jet Propulsion Laboratory DE-200 
[Standish, 1982]. 

A tabulation of the angle and the distance of the planets from the Earth is 
provided in Table 6.1, as given by the Astronomical Almanac for the January 1, 1984, 
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epoch date. This table provides the initial locations of the planets in a heliocentric 
coordinate system which enables an approximate determination of the relative location 
to the Earth. Clearly, the closer a planet is to the Earth, the larger its direct effect on the 
satellite’s motion can be, but the actual magnitude is also dependent on the planets 
mass. 


Table 6.1 

Heliocentric coordinates for the planets 
January 1, 1984 


Planet 

longitude (°) 

latitude (°) 

distance from sun 

Mercury 

97.05 

5.17 

0.310727 

Venus 

178.41 

3.19 

0.719986 

Earth 

302.55 

18.94 

0.98400 

Mars 

169.45 

1.35 

1.66087 

Jupiter 

263.22 

0.22 

5.28589 

Saturn 

219.12 

2.23 

9.83198 


6.6 Luni-Solar Effects 

The initial configuration of the Sun and the Moon relative to the orbit plane of 
the GRM satellites will be a significant factor and will influence the overall magnitude 
of the range-rate. The forces due to the Sun and the Moon are also described by 
Equation (6.1), where the Sun and Moon are taken as the perturbing bodies. The 
effects of the Sun and the Moon on the relative range-rate were investigated separately 
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and also combined. 

The largest, single temporal perturbation on the satellites relative range-rate is 
due to the Moon (Fig 6.6a). In the formulation used, the Moon is assumed to be a 
point mass. Since the simulation to determine the Moon’s effect is 32 days, the 
satellites are exposed to an entire revolution of the Moon about the Earth. The 
maximum perturbation encountered in this simulation occurred when the Moon is near 
the inertial Y-axis. 

The solar effects on the relative range-rate are presented in Figure 6.6b. Like 
the Moon, the Sun has its greatest affect when the Sun is along the inertial Y-axis. 
Since the mission is to last six months, this configuration will be encountered at least 
once, and perhaps twice. For the January 1 epoch date, the Sun is initially near the 
longitude of 279°, which is very close to the negative Y-axis. It is possible to decrease 
the effect the Sun will have on the orbital motion of the satellites by permitting the Sun 
to be on the Y-axis only once. This can be achieved by beginning the mission in March 
or September. 

The effects of the Moon and the Sun are combined to be the luni-solar effect 
(Figure 6.6c). In a paper by Estes and Lancaster [1976b], it was stated that the luni- 
solar effects can be minimized by placing the satellites in a orbit plane that is 
perpendicular to both the ecliptic plane and to the equatorial plane. The orbit plane that 
satisfies both of these criteria is the plane that is perpendicular to the vernal equinox (the 
X-axis) which is the Y-Z plane. The GRM satellites are initially in the Y-Z plane, 
therefore, in the plane of minimum luni-solar effect. This plane should precess at the 
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rate the vernal equinox precesses (360°/26000 years), but since the mission is only six 
months long and this rate is small and can be neglected. 

The perturbations of the Moon's effect on the Earth's oblateness were included 
in this section. The force equation is provided by Equation (6.1), except that the 
indirect term in the force equation, t-Jt?, is replaced by the indirect effect due to 
oblateness, VU(r;)/(GM) [Moyer, 1971] where VU(r) is the gradient of the Earth 
geopotential field with respect to rj, the distant from the Moon to the Earth. The relative 
range-rate plot (Figure 6.6d) contains a periodic effect whose magnitude remains under 
the ±1 pm/s level. Since there appears to be no secular growth in the amplitude over 
time, this effect on the satellites' motion can be excluded from the dynamic model. 

6.7 Relativistic Effects 

The model used in UTOPIA to calculate the relativistic effects on a satellite's 
motion assumes that the spacecraft is a massless particle revolving around a point mass 
[Moyer, 1971]. The dominant effect of relativity is the effect on the motion of perigee. 
Because of the low altitude of the satellites in this mission, the relativistic effects on the 
satellites’ motion were expected to be significant. Using mean orbital elements, the 
perigee advance rate is approximately 0.0633°/day. According to general relativity, the 
contribution of the relativistic perturbation to the equations of motion is: 


r = 4p/(c 2 r 3 ) { [ p/r - r • r ] r + ( r • r ) r } 


where c is the speed of light in a vacuum, p is the gravity parameter of the Earth, r is 
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the position vector, and f is the velocity vector of the satellite [Moyer, 1968]. 

Figure 6.7 illustrates that the relativistic effects influence the satellites' relative 
range-rate above the ±1 |im/s level. The amplitude of the relativistic effects increases to 
±16 pm/s in 32 sidereal days, indicating that the Newtonian model is insufficient to 
model the acceleration of these satellites. 

6.8 Effects of the Perturbations on the Initial Conditions 

The results of the combination of all the major temporal effects on the relative 
range-rate are illustrated in Figure 6.8a; the Moon’s effect on the Earth oblateness and 
all the planets, except for Jupiter, were excluded from the model. The relative motion 
plots are provided in Figures 6.8b and 6.8c. The relative range (Figure 6.8b) indicates 
that a drift between the satellites of 10.95 m/day has been incurred due to the additional 
forces. The satellites' groundtrack repeat has also been affected; primarily, the 
longitudes of the two satellites was west of the 32 sidereal day closure point. For the 
nominal case, the two satellites had a groundtrack closure of within 100 meters; with 
the temporal perturbations, the satellites close to within four kilometers. 

A correction to the nominal initial conditions was made using the technique 
described in Chapter 4. However, this is a linear technique and is not completely suited 
for these temporal perturbations. From Table 4.1, if there is a total error in latitude of 
0.002°, then the drift rate should be approximately 3.5 m/day. Also, the error in 
longitude is predicted to be negligible. This is not the case with the temporal effects. 
Instead, with the same error in latitude the drift rate was almost three times greater, and 
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the error in longitude was close to 0.2° for each satellite. The correction for the initial 
conditions that will result in a more accurate groundtrack closure had to be obtained 
iteratively, that is, corrections to the initial conditions were determined and the resulting 
closures and drift rate were calculated. If these conditions were not acceptable, the 
procedure was performed again. Nutation is the cause of the nonlinearity, since the 
coordinate system will nutate nearly two degrees in 32 sidereal days. 

After two iterations, the final values for latitude were equal to the initial values. 
The corrections to the initial conditions were 0.3544676 m for the leading satellite and 
0.31 17409 m for the trailing satellite, both in the negative z direction. Longitude was 
still west of its desired value but closer to an exact repeat. The drift rate was reduced to 
1.44 m/day and the closure was under three kilometers due to the remaining longitude 
error. The relative motion is illustrated in Figures 6.9a and 6.9b. In addition, the 
combined effect of all the temporal perturbation with the new initial conditions was 
determined (Figure 6.9c). There was no discernible improvement in the relative range- 
rate; that is, changing the initial condition slightly does not seem to significantly alter 
the perturbation effects on the satellites' relative motion. 

6.9 Effect of the Perturbations on the Frozen Orbit 

An investigation of the perturbation effects on the character of the frozen orbit 
was considered. Figures 6.10a and Figure 6.10b indicate that the perturbed orbits will 
remain frozen since the phase plane diagram and the eccentricity versus the argument of 
perigee were within the same patterns of the unperturbed 9x9 geopotential field frozen 
orbit (Figures 3.4a and 3.4b). The frozen orbit characteristics were not significantly 
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influenced by the temporal perturbations investigated in this study. 

The luni-solar effects on the frozen orbit were investigated for the proposed, 
Navy satellite NROSS [Cefola, et al., 1986] and for SEAS AT [Nickerson, et al., 
1978]. Their results indicated that the luni-solar effects do not alter the frozen orbit, at 
least for the time period considered. These studies were made on satellites with much 
higher altitudes than is planned for GRM, indicating that the luni-solar effects should 
not interfere with the GRM frozen orbit characteristics as was illustrated in Figure 6.10. 

6.10 Summary 

With the accuracy levels required of this mission, the dynamical model will 
need to be detailed and complete. Any perturbations to the satellites’ orbits that could 
alter the relative range-rate will have to be accounted for in the modeling in order to 
correctly identify the geopotential field. In some cases, for instance ocean tides, an 
error in the model of only 10% may produce signals exceeding the measurement 
precision which would be detrimental to the recovery the the Earth's gravity field. 
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e is the obliquity of the ecliptic 
Z is the pole of the ecliptic 
z is the rotation axis of the Earth 


Figure 6. 1 

Precessional and nutational motion 
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DIFFERENCE IN RELATIVE RANGE-RATE 
OCEAN TIDAL EFFECTS 
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Figure 6.4 
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DIFFERENCE IN RELATIVE RANGE-RATE 
ALL TEMPORAL EFFECTS INCLUDED 
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Figure 6.8a 


Difference in relative range -rate 
All temporal effects combined 


36. o 


185 


INITIAL CONDITIONS FOR 9X9 GEOPOTENTIAL 
ALL TEMPORAL PERTURBATION INCLUDED 
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Figure 6.8b 
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Figure 6.8c 
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ADJUSTED INITIAL CONDITIONS FOR 9 X 9 GEOPOTENTIAL 
ALL TEMPORAL PERTURBATION INCLUDED 
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Figure 6.9a 

Relative range: adjusted initial conditions 
All temporal effects combined 
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Figure 6.9b 

Relative range-rate: adjusted initial conditions 
All temporal effects combined 
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DIFFERENCE IN RELATIVE RANGE-RATE 
ADJUSTED: ALL TEMPORAL EFFECTS INCLUDED 
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Figure 6.9c 

Difference in relative range-rate 
All temporal effects combined 
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Figure 6.10b 

Eccentricity versus argument of perigee: Short periodic effects were included. 
Mean orbital elements were frozen. Temporal effects included 





CHAPTER 7 


CONCLUSIONS AND SUMMARY 

This study concentrated on the pre-mission and post-mission phases of the 
proposed Geopotential Research Mission. The pre-mission phase determined a strategy 
for calculating a set of initial conditions which required an entire repeat cycle that met 
the mission specifications of a frozen, polar orbit with repeating groundtracks for both 
satellites. Corrections to the initial conditions were determined using the same strategy 
after only one week of mission time with little loss of accuracy. The post-mission 
phase determined a nominal set of initial conditions along with a reduced geopotential 
field to produce an orbit that satisfied the mission requirements for orbit accuracy. 

7.1 Summary 

The definition and usefulness of a frozen orbit was discussed in Chapter 3. For 
a polar, frozen orbit, the mean argument of perigee location is a constant. It was 
demonstrated that frozen orbits can maintain a repeating groundtrack more easily than a 
non-frozen orbit. Once obtained, the frozen orbit configuration is very stable, and 

i 

1 perturbations as large as 1000 meters in the orbit position did not destroy the integrity 

of the orbit's characteristics. 
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A strategy for determining a set of initial conditions was described in Chapter 4. 
Once a set of initial conditions met the criteria for a frozen orbit, an adjustment based on 
a linear calculation to these initial conditions will allow them to have a repeating 
groundtrack after the specified number of days. The adjustments needed in the initial 
conditions were independent of the geopotential field that influenced the satellites' 
motion. Two methods were provided, one that required the entire repeat cycle (in this 
study, 32 sidereal days) to determine the proper adjustments, and another that used a 
maximum of only one week of mission time. The sensitivity of the initial conditions to 
orbit insertion errors and the range of linear reliability was also investigated. 

A simulation of the satellites was performed and the results of this simulation 
were discussed in Chapter 5. The simulation used an Ohio State University 
geopotential field, which consisted of a 180 x 180 field plus coefficients to degree 300 
and up to order 10 [Rapp, 1981]. The simulation showed that the satellites' relative 
motion was highly influenced by certain resonant terms, particularly orders 82, 33, 49, 
164, 16, and 17. The ephermerides of each satellite were used as a set of observations 
to simulate the orbit determination process. To reduce the difference between the 
observations and the nominal trajectory based on the GEM10B geopotential field, the 
first two pairs of each of the resonant coefficients plus the zonal harmonics J2 and J3 
were estimated. With this nominal geopotential field, the nominal orbit accuracy was 
reduced to satisfy the gravity mission specifications. 

Chapter 6 provided a study of perturbation effects. The effects chosen for this 
investigation were: precession, nutation, and polar motion, planetary, luni-solar, 
relativity, solid Earth tides, and ocean tides. Except for polar motion, these 
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perturbations were shown to influence the satellites' relative range rate in excess of the 
±1 jim/s requirement, and therefore, must be accounted for in the modeling process. 

This study has been based on an epoch date of January 1, 1984. A different 
initial epoch or different mission requirements, such as groundtrack repeat frequency or 
separation distance, will require a different set of initial conditions. However, the 
procedure provided in Chapter 4 to obtain the initial conditions remains valid regardless 
of the final mission requirements. Some of the results, such as in the temporal effects 
on the satellites' relative motion (Chapter 6), will depend upon the epoch date selected. 
In addition to the epoch date, the results from the simulation are dependent upon the 
geopotential field used to generate the satellites' trajectories. Another simulation will 
require a different nominal trajectory that will meet the mission residual requirements. 
Also, additional harmonic terms may need to be included in the nominal model. 
However, the relative power of the resonance terms should remain as presented for the 
same satellite altitude of 160 km that was selected for this study. 

7.2 Additional Research 

Additional studies proposed for the Geopotential Research Mission include the 
following: 

o Expanding the geopotential field to a full 360 x 360 field for the simulation. 
Some of this work has already been done and the results of this larger 
simulation were partially discussed in Chapter 5. 
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o A study of the atmospheric effects on the outer shell, and the manner in which 
thrusts associated with the compensations for the nongravitational force 
influences the satellites' relative motion. 

o An investigation of tracking systems, such as TRANET and GPS and their 
ability to provide a sufficiendy precise determination of the orbits. 

o The investigation of techniques for the recovery of the coefficients of the 
geopotential field with simulated data. 

o A simulation of only one satellite equipped with a gravity gradiometer instead of 
the dual satellite configuration to be used to recover the geopotential field. 


APPENDIX A 


Lagrange's Planetary equations for the modified set of orbit elements described 
in Chapter 3 are as follows [Taff, 1978]: 

d = 2 (a /p ) 1/2 dR/da 

fl = (\-e 2 )/L dR/d£ - fyoti /L dR/di - i](l-e 2 )/L[l+(l-e 2 ) 1 / 2 ]dR/da 
4 = -(l-e 2 )/L dR/dr| + T] coti /L dR/di - £(l-c 2 )/L[l+(l-<? 2 ) 1 / 2 ]dR/dc 
di/dt = L'^coti £ dR/dri - r\ dR/d£ + dR/da) - cos i dR/dQ] 

Q = esc i IL dR/d/ 

d = (l-e 2 )/L [l+(l-e 2 ) ly2 ] (iidR/dn+^dR/d^) - 2(a /p) 1 / 2 dR/da - coti /L dR/d i 
where 

L 2 = pa ( 1 - e 2 ) 

Disturbing function, R, used includes J 2 , J 3 , J 5 , J 7 , and J9 only. 

R = V SP + Vjjp +V$ec 
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Short period contribution to the disturbing function: 

V S p = -3 pa e 2 /a J 2 [sin 2 /(7/8^cosa + 5/8risina - 7/8(^cos3a + 
Tjsin3a) - l/4cos2a) - l/2(^cosa + Tjsina)] 

where a = M + to 

a = to + ntp 


Long period contribution to the disturbing function [Kaula, 1966]: 

V LP = -3/2 (p/a 3 ) 172 (ae/a) 3 T] sin i {J 3 (l-5/4sin 2 / ) - 5/2 J 5 (a Ja) 1 ( 1 7/2sin 2 /+2 l/8sin 4 /) 
+ 35/8 J 7 (a e/a) 4 (l-27/4sin 2 / + 99/8sin 4 / - 429/64sin 6 /) - 105/16 J 9 (a <Jaf (1 - 
11 sin 2 / + 143/4sin 4 / - 715/16sin 6 / + 2431/128sin 8 /)} 

Secular contributions to the disturbing function: 

V S EC = -h V a e 3 /a 3 (3/4sin 2 / - 1/2) / (l-e 2 f n 

Cook's analysis excluded the short period contributions to the disturbing 
function, V SP . By excluding the short period term. Cook was able to find a analytic 
solution to the equations of motion for £ and f|. For Cook's solution 3R/9a is zero. 


I 


£ = -kri + C 
f|=k£ 
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From Cook's [1966] solution: 

£ = A cos ( kt + a ) 

T| = A sin ( kt + a ) + C/k 

where 
C = Vlp/ti 

k = 3 ( n/a 3 ) 1/2 J 2 (ae/a) 2 (l-5/4sin 2 /) - 

5 Oi /a 3 ) m J 4 (a e/a) 4 (105/64 sin M - 15/8 sin 2 / + 3/8) 


Note that: 

e = (Z ) 2 + t| 2 ) 1/2 

co = tan'^Tj / £) 



APPENDIX B 


Program FIXDRF calculates the correction needed in each of the satellite's 
initial conditions that will eliminate the drift between the satellites and that will insure a 
closure of within one meter. 

PROGRAM FIXDRF ( INPUT, OUTPUT, TAPE5=INPUT, TAPE6=OUTPUT ) 

C 

C THIS PROGRAM FIXES THE DRIFT BETWEEN TWO SATELLITES GIVEN THE 

C GRAVITY PARAMETER GM, AND THE MEAN EARTH RADIUS AE . FIRST ORDER 

C ASSUMPTION ONLY, NO GRAVITY COEFFIECIENTS ARE INCLUDED. THE ORBIT 

C HAS TO BE NEARLY CIRCULAR AND POLAR. ADJUSTMENT IS IN THE Z POSITION. 

C 

REAL NBAR, MU, J2, NRATE 

COMMON / CHANGE / DLONGl , DLONG2 , DLAT1, DLAT2 
C 

C MEAN ORBITAL ELEMENTS 

A = 6523600.233433 
E * .00153496544 
J2 = .00108262808458 
AE = 6378137. 

MU = 3 . 986004 4E14 

NBAR = SQRT ( MU / A**3 ) 

NRATE = (A * ( 1.- E**2)/AE ) **2 * 2./3./ < NBAR * J2 ) 

C 

C LATITUDES OF EACH SATELLITE INITIAL AND FINAL 

PHI 1 = 88.688 
PHI1F * 87.3059 
PHI2 = 88.69072 
PHI2F = 89.87647 
C 

C IF DLAT IS GREATER THAN 0 SAT IS BEHIND 

C IF DLAT IS LESS THAN 0 SAT IS AHEAD 

C CHANGE IN LATITUDES 

DLAT1 = PHI1F - PHI 1 
DLAT1 - .015 
DLAT2 = PHI2 - PHI2F 
DLAT2 = .003 
C 

WRITE (6,220) 

WRITE (6,200) DLAT1 , DLAT2 
C 

C RADIUS 

R = A 
C 

C DRIFT RATE IN METERS PER DAY 

DR = 9.6469 
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C CHANGE IN MEAN MOTION IN RADIANS PER SECOND 

DN « DR / (R * 86400.) 

C 

WRITE (6, 140) 

WRITE (6,100) DN 

C CHANGE IN Z POSITION IN METERS DUE TO DRIFT 

DZ = DN * ( 2./3. ) * A / NBAR 
WRITE (6, 150) 

WRITE ( 6, 100) DZ 

C CALCULATE CHANGE IN TWO SATELIITES Z POSITIONS TO ADJUST LATITUDE 

C 

C CHANGE IN LATITUDE POSSIBLE DUE TO DRIFT RATE 

DLAT = DN * 2757250.896 * 57.29577951 
WRITE (6, 100) DLAT 
C 

CALL DRFADJ ( DLAT1, DLAT2 , DLAT, DZ, DZ1, DZ2 ) 

C 

WRITE (6,250) 

WRITE (6, 200) DZ1 , DZ2 
C 

STOP 

100 FORMAT ( 5X, E14 . 7 ) 

140 FORMAT ( 5X, 1 CHANGE IN MEAN MOTION ' ) 

150 FORMAT ( 5X, * CHANGE IN Z DUE TO DRIFT' ) 

200 FORMAT ( 5X, 2E14.7) 

220 FORMAT ( 5X, ' CHANGE IN LATITUDE ' ) 

230 FORMAT ( 5X, 1 CHANGE IN LONGITUDE » ) 

250 FORMAT ( 5X, ' CHANGE IN Z FOR EACH SATELLITE ' ) 

300 FORMAT (5X, 6E2 1 . 13 , / , 5X, 6E21.13) 

END 

C 

SUBROUTINE DRFADJ ( DLAT1 , DLAT2 , DLAT, DZ, DZ1, DZ2 ) 

KFLAG = 1 
Q = -.5 

IF ( DLAT1 .GE. 0. ) Q = .5 
DZ1 - 0. 

DZ2 - 0. 

IF (ABS (DLAT2) ,GT . ABS (DLAT1 ) ) KFLAG * -1 
IF { KFLAG .EQ. 1 ) GO TO 10 
DLAT2 = ABS (DLAT2) - ABS (DLAT) 

DLAT2 = ABS (DLAT2 ) 

DZ2 = ABS (DZ) 

GO TO 20 
C 

10 DLAT1 = ABS (DLAT1 ) - ABS (DLAT) 

DLAT1 = ABS (DLAT1 ) 

DZ1 = ABS (DZ) 

20 CONTINUE 

DELAT = DLAT2 

IF ( DLAT 2 .GT. DLAT1 ) DELAT = DLAT1 
DZ1 « Q * ( DZ * DELAT / DLAT + DZ1 ) 

DZ2 = Q * ( DZ * DELAT / DLAT + DZ2 ) 

RETURN 

END 
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